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1  Introduction 

The  assignment,  or  matching,  problems,  such  as  the  Linear  Assignment  Problem  (LAP)  and  the 
Quadratic  Assignment  Problem  (QAP),  are  of  theoretical  and  practical  importance  in  a  number  of  fields, 
from  discrete  mathematics  and  computer  science  to  archeology  and  sports  (for  details  on  theory  and 
applications  of  assignment  problems  see,  among  others,  Pardalos  and  Pitsoulis,  2000;  Pentico,  2007; 
Burkard  et  ah,  2009).  The  underlying  combinatorial  structure  of  the  feasible  set  common  to  these  prob¬ 
lems  is  that  of  a  perfect  matching  on  a  bi-partite  graph,  i.e.,  a  bijective  mapping  between  the  elements  of 
two  sets.  The  present  endeavor  is  concerned  with  generalizations  of  the  classical  assignment  problems 
where  the  underlying  combinatorial  structure  is  based  on  hypergraph  matchings. 

A  hypergraph  is  a  pair  %  =  (V,  £)  where  V  is  the  set  of  vertices,  and  £  is  the  set  of  hyperedges; 
a  hyperedge  e  €  £  is  a  subset  of  V  containing  two  or  more  vertices  (Berge,  1989;  Bollobas,  1998). 
Similarly  to  regular  graphs,  a  hypergraph  H  is  called  /c -partite  if  its  vertices  can  be  partitioned  into  k 
independent  subsets  such  that  no  vertices  within  a  subset  are  connected  by  an  edge.  If  all  edges  of  a 
hypergraph  contain  exactly  r  vertices,  the  hypergraph  is  called  r -uniform. 

In  this  paper  we  restrict  our  discussion  to  d -partite  d -uniform  hypergraphs  'Hj\n  =  (V,  £),  whose  sets 
of  vertices  V  arc  partitioned  into  d  independent  sets  V/,  each  of  cardinality  n:  V  =  Uy=i  V/,  |  V/ 1  =  n. 
Consequently,  each  hyperedge  of  Hd\n  contains  exactly  one  vertex  from  each  independent  set  V/,  and 

thus  it  can  be  represented  as  a  vector  (z  i . ij)  e  Vi  x  ■  •  •  x  Vj.  This  convention  allows  us  to  denote 

the  vertices  of  each  Vjt  as  1 , 2, so  that  a  hyperedge  can  be  presented  as  (z  1 . i^)  €  {1 . n}d . 

Finally,  it  is  assumed  that  the  d -partite  d -uniform  hypergraph  'Hd\n  is  complete,  i.e.,  it  contains  nd 
hyperedges. 

Then,  a  (perfect)  matching  \i  on  'Hci\„  is  formed  by  a  set  of  n  hyperedges  that  do  not  share  any  vertices, 
or,  equivalently,  which  have  the  property  that  each  vertex  of  the  hypergraph  belongs  to  exactly  one 
hyperedge: 

p  =  {{e\, . . . ,  e„j  |  et  €  £ ,  <?,  n  q  =  0,  1  <  z,  j  <  zj}. 
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Since  a  hyperedge  of  %d\n  contains  exactly  one  vertex  from  each  of  d  independent  sets  V/c,  it  can  be 

represented  as  a  vector  {i\, . ic[ )  e  [«]  ,  where  the  set  [[«]  =  {1 . n}  is  used  to  label  the  vertices 

of  each  V/t-  Then,  the  set  j\4('H(j\n)  of  perfect  matchings  on  'H(j\n  can  be  represented  in  mathematical 
programming  form  as 

M{Ud\n)  =  ^xe{0,l}nd 

Are'i</]\{r} 


Y  xh...id  =  \,  ir  e  In},  r  g  [rf]  j,  (1) 

4e[n]  ) 


where  vy  , ...  id  —  1  if  the  hyperedge  (i\ . i(j)  is  included  in  the  matching,  and  x, , ...  ,d  =  0  otherwise. 

If  the  cost  of  hypergraph  matching  fi  is  given  by  T>(/i),  the  general  combinatorial  optimization  problem 
on  hypergraph  matchings  can  be  formulated  as 


min  <  0  ( /x) 


li  e  M{Hd\n)  |  > 


(2) 


where  M  is  the  set  of  all  perfect  matchings  on  the  hypergraph  Hci\n-  A  mathematical  program¬ 

ming  formulation  of  the  general  problem  (2)  can  be  written  as 


min  <3>(x) 

n 

n 

(3a) 

s.t.  £•• 
*2=1 

n 

•  E  xi\  —  id  =  1> 

id  =  1 

n  n  n 

h 

=  1,.. 

.  ,n, 

(3b) 

E  •• 

*i=t 

n 

•  E  E  E  xi\-id  = 

*Ar— 1  =1  4-+l  =  l  4  =  1 

n 

ik 

=  1,.. 

. ,  n,  k  —  2,. 

. . ,  d  —  1, 

(3c) 

E  •• 

*i=t 

•  E  xi\  —  id  —  '  • 

4-i  =  l 

id 

=  1, .. 

. ,  n , 

(3d) 

X  e  {0,l}nrf, 

(3e) 

where  xq ...  id  —  1  if  the  hyperedge  (/ ] . i(j)  is  included  in  the  matching,  and  xq ...  l(,  =  0  otherwise. 

The  mathematical  programming  formulation  (3)  of  the  hypergraph  matching  problem  (2)  makes  it  natural 
to  call  problem  (3)  a  multidimensional  assignment  problem  (MAP),  where  d  stands  for  the  number  of 
“dimensions”,  and  n  is  the  number  of  “elements  per  dimension”  (recall  that  d  is  equal  to  the  number  of 
independent  subsets  of  vertices  in  the  hypergraph  and  n  denotes  the  number  of  vertices  in  each 
Vk)-  In  the  sequel,  we  will  use  the  terms  “combinatorial  optimization  problem  on  hypergraph  matchings” 
(or  “hypergraph  matching  problem”  for  short)  and  “multidimensional  assignment  problem”  in  reference 
to  (2) — (3)  interchangeably. 

Depending  on  the  particular  form  of  d>(-),  a  number  of  combinatorial  optimization  problems  on  hyper¬ 
graph  matchings  can  be  formulated.  In  this  paper,  we  analyze  several  problems  that  arise  as  special  cases 
of  (2)-(3)  when  the  matching  cost  function  <f>  has  the  form 

=  Li  (4) 

9 •  •  •  ?  eim 

where  JJ  is  an  operator  defined  over  some  set  of  cost  elements  {0}  indexed  by  the  hyperedges  of 
For  instance,  if  the  cost  function  in  (2)  is  defined  as  d>(/z)  =  (j>e ,  or,  equivalently,  as  a  linear  form  over 
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variables  xIr..Jrf  in  (3),  one  obtains  the  so-called  Linear  Multidimensional,  or  Multi-index  Assignment 
Problem  (LMAP):1 


Z 


d\n  — 


min 
s.  t. 


n  n 

H  '  '  '  H  I  -  id  Xl  I  -  id 
*1  =  1  id  =  1 


(3b)-(3e). 


(5) 


Taking  <I>(x)  as  a  quadratic  form  over  x  e  {0,  1  \n'' ,  or,  in  other  words,  d>(/x) 
the  Quadratic  Multidimensional  Assignment  Problem  (QMAP): 


L  (Pe,ej,  leads  to 
e, ,  ej  e/i 


Qd\n 


min 
s.  t. 


n  n  n  n 

H  ' ' '  H  H  H  Qii-idji-jd  xh-id  xh-id 

*1  =  1  id  =  l7i  =  l  id  =  1 


(3b)-(3e). 


(6) 


Similarly,  if  one  is  interested  in  minimizing  the  largest  element  of  the  corresponding  linear  or  quadratic 
forms,  the  Linear  Bottleneck  MAP 


Wd\n  = 


min  max  d>iv..id  x, 

*l,  —  .*rf6{l,  — .«} 

s.  t.  (3b)-(3e) 


(V) 


■  .  max  .  Qii-idji-jd  xii-id  Al\-Jd 

*1, •••*</ 


Xu 


(8) 


and  Quadratic  Bottleneck  MAP 
Ud\n  =  min 

s.  t.  (3b)-(3e) 

problems  are  obtained. 

A  matching  /x  =  { ( f { 1  ^ . z^), . . . ,  (z'f^ . z^)}  on  %d\n  can  be  conveniently  presented  in  the 

matrix  form. 


/  i\V)  i2 


H  = 


(i) 


;(2)  -(2) 


id)  \ 

f2) 


1 


V 


:(«)  :(«) 


!  1 


;(n) 


(9) 


/ 


where  each  column  (ill  if2) 


:  (**)  \  T 


,  ik  i'k)  ,  k  =  1,  •  •  • ,  d,  is  a  permutation  of  the  set  [[«].  Since  the  rows 

of  matrix  (9)  represent  hyperedges  of  'Hj\n.  they  can  be  permuted  so  that  in  the  /< -tli  dimension  we  have 


.  h 
of  (2): 


)  =  (1 . zz ).  This  allows  for  a  permutation  representation  of  a  feasible  solution 


P  = 


/ 1 

7Ti(l)  •• 

•  *d- t(l)  ^ 

2 

TTl  (2)  ■■ 

•  *d- 1  (2) 

V  n 

ni(n)  •• 

•  nd-i(n)  / 

(i,  JTl,  .  .  .  ,  7 X d  —  1), 


(10) 


1  Problem  (5)  is  also  known  in  the  literature  as  the  axial  MAP  and  is  distinguished  from  the  so-called  planar  MAP,  whose 
underlying  combinatorial  structure  is  based  on  latin  squares;  see,  e.g.,  Burkard  et  al.  (2009). 
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where  each  :  [«]  i->  [«]  is  a  permutation  of  the  set  [«],  and  i  is  the  identical  permutation,  i(i  )  =  i, 
i  =  1 . n. 

The  multidimensional  assignment  problems  (5)— (8)  also  admit  useful  permutation  formulations;  for  in¬ 
stance,  the  Lineal-  MAP  (5)  can  equivalently  be  written  as 

n 

Zd\n  —  min  yi  3T1  Q~)— JTrf— 1 0')’ 

7Tl  ,...,7?^/— 1  €11/7  , 

1  =  1 

where  is  a  permutation  of  the  set  [[«],  and  Y\n  is  the  set  of  all  such  permutations.  Problems  (6)— (8) 
can  be  reformulated  analogously.  The  permutation  formulation  makes  it  easy  to  see  that  the  cardinality 
of  the  feasible  set  of  (2)  or,  equivalently,  (3),  is 

\M(J~Ld\n)\  =  n\d~l. 

The  Lineal'  MAP  (5)  was  first  introduced  by  Pierskalla  (1968),  and  had  found  numerous  applications  in 
the  areas  of  data  association,  sensor  fusion,  multi-sensor  multi-target  tracking  (Murphey  et  al.,  1998a; 
Poore,  1994a;  Kirubarajan  et  ah,  2001;  Poore  and  Gadaleta,  2006;  Pusztaszeri  et  ah,  1996;  Andrijich 
and  Caccetta,  2001;  Chummun  et  al.,  2001),  image  recognition  (Veenman  et  al.,  1998),  and,  recently, 
peer-to-peer  refueling  of  space  satellites  (Dutta  and  Tsiotras,  2008);  for  a  detailed  discussion  of  the 
applications  of  the  MAP,  see,  e.g.,  Burkard  and  Cel  a  (1999a)  and  Burkard  (2002).  A  three-dimensional 
version  ( d  —  3)  of  the  Quadratic  MAP  (6)  was  considered  by  Sarni  a  et  al.  (2005)  and  Hahn  et  al.  (2008) 
in  application  to  design  of  robust  wireless  transmission  systems. 

The  objective  of  the  present  paper  was  to  elucidate  the  properties  of  large-scale  combinatorial  optimiza¬ 
tion  problems  on  hypergraph  matchings  (5)— (8)  by  analyzing  the  random  instances  of  the  corresponding 
problems,  i.e.,  by  assuming  that  the  assignment  costs  ld...  in  (5)— (8)  are  iid  random  variables  from 
a  specified  probability  distribution.  The  introduced  multidimensional  assignment  problems  (5)— (8)  con¬ 
stitute  the  primary  focus  of  the  present  paper;  as  it  will  be  seen  below,  the  proposed  techniques  admit 
extension  to  other  forms  of  cost  function  <J>  given  by  (35)  in  the  general  formulation  (2) — (3). 


Related  work  Studies  of  random  instances  of  combinatorial  problems  based  on  bipartite  graph  match¬ 
ings  that  are  given  as  special  cases  of  (2)— (3)  with  cl  —  2  have  been  an  active  area  of  research  for  the  last 
several  decades.  A  number  of  important  results  have  been  obtained  in  the  scope  of  random  linear  and 
quadratic  assignment  problems,  some  of  them  within  the  last  few  years;  for  a  detailed  review  of  random 
assignment  problems  see  Krokhmal  and  Pardalos  (2009). 

Considerable  attention  has  been  paid  to  studies  of  the  limiting  behavior  of  the  expected  optimal  value  of 
the  Lineal'  Assignment  Problem,  a  special  case  of  (5)  with  cl  =  2.  In  the  works  of  Walkup  (1979);  Karp 
(1987);  Lazarus  (1993);  Goemans  and  Kodialam  (1993);  Olin  (1992);  Coppersmith  and  Sorkin  (1999), 
a  series  of  successively  tighter  upper  and  lower  bounds  for  the  expected  cost  E[Z2\n]  of  LAP  have  been 
obtained  in  the  assumption  that  the  assignment  costs  are  iid  uniform  on  [0,  1]  or  exponential  with  mean 
1.  Bounds  on  E[Z2\n]  for  general  cost  distributions  have  been  furnished  in  Frenk  et  al.  (1987),  and  Olin 
(1992)  presents  bounds  for  Z2 \n  that  hold  with  probability  1. 

Perhaps,  the  most  widely  known  results  concerning  the  random  LAP  are  the  conjectures  by  Mezard  and 
Pai'isi  (1985)  and  Parisi  (1998)  that  the  expected  optimal  cost  of  a  random  LAP  with  iid  uniform  [0, 1]  or 
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exponential  with  mean  1  assignment  costs  satisfies 


lim  E[Z2i„]  =  —  1.645  and,  moreover,  E[Z2iJ  =  (12) 

n  — ►  oo  6  z— ^  k 1 

k=  1 

The  existence  of  the  limit  in  (12)  was  argued  by  Mezard  and  Parisi  (1985)  using  a  non-rigorous  replica 
method  and  was  experimentally  observed  in  Donath  (1969)  and  Pardalos  and  Ramakrishnan  (1993),  but  a 
theoretical  proof  of  the  limit  in  (12)  was  first  furnished  by  Aldous  (1992,  2001).  Within  the  last  five  years, 
the  second  conjecture  in  (12)  for  finite-sized  LAPs  due  to  Parisi  (1998)  has  been  proven  independently  by 
Linusson  and  Wastlund  (2004)  and  Nair  et  al.  (2005),  and  its  generalizations  (Coppersmith  and  Sorkin, 
1999;  Linusson  and  Wastlund,  2000;  Buck  et  ah,  2002)  for  the  case  of  random  LAPs  with  non-square 
matrices  (the  so-called  Ar-assignment  problem)  have  been  proven  in  Wastlund  (2005a,b). 

The  methods  of  probabilistic  analysis  have  played  a  pivotal  role  in  explicating  the  computational  proper¬ 
ties  of  many  hard  combinatorial  optimization  problems,  and,  particularly',  the  Quadratic  Assignment 
Problem  that  is  given  by  (6)  with  d  =  2.  While  the  QAP  is  known  to  be  NP-complete  and  non- 
approximable  (Sahni  and  Gonzales,  1976)  and  thus  very  difficult  to  solve  exactly  (currently,  instances 
of  the  QAP  of  sizes  up  to  n  =  30  can  be  solved  rotinely,  see  Anstreicher,  2003;  Loiola  et  al.,  2007),  it 
was  noticed  that  good  quality  sub-optimal  solutions  of  the  QAP  can  be  obtained  relatively  easily  with 
many  heuristic  algorithms.  It  turns  out  that  this  strange  behavior  can  be  explained  within  the  probabilis¬ 
tic  framework,  i.e.,  by  analyzing  the  random  instances  of  the  QAP.  In  particular,  Burkard  and  Fincke 
(1982a)  were  the  first  to  point  out  that  in  a  random  QAP  the  ratio  between  the  minimum  (optimal)  cost 
02\n  and  the  maximum  cost  Q2<1’1  approaches  1  in  probability  as  the  size  of  the  problem  increases: 

lim  pj-^ —  <  1  +  0(n~0  225)  |  =  1.  (13) 

n-+oo  (  Q2 1„  ) 

Frenk  et  al.  (1985)  and  Rhee  (1988,  1991)  established  stronger  versions  of  (13),  including  a  proof  that 
the  cost  ratio  in  (13)  converges  to  unity  not  only  in  probability,  but  almost  surely  (a.s.).  It  must  be  noted 
that  the  re  mark  able  observation  of  Burkard  and  Fincke  (1982a)  has  led  to  the  discovery  of  an  entire  class 
of  combinatorial  optimization  problems  with  similar  behavior  (Burkard  and  Fincke,  1985;  Szpankowski, 
1995;  Albrecher  et  ah,  2006). 

A  number  of  results  along  the  lines  of  (12)  and  (13)  have  been  obtained  for  other  combinatorial  optimiza¬ 
tion  problems  based  on  bipartite  graph  matchings,  such  as  linear  and  quadratic  bottleneck  assignment 
problems,  bi-quadratic  assignment  problem,  etc.,  in  Pferschy  (1996);  Burkard  et  al.  (1994);  Burkard  and 
Fincke  (1982b);  Albrecher  (2005),  and  others;  see  Krokhmal  and  Pardalos  (2009)  for  details. 

As  regards  the  combinatorial  optimization  problems  on  hypergraph  matchings  (2)-(3),  relatively  few 
attempts  have  been  made  in  the  literature  to  investigate  their  properties  in  large-scale  instances.  Dyer, 
Frieze,  and  McDiarmid  (1986)  derived  an  upper  bound  on  the  expected  optimal  cost  E[Z^f  ]  of  the  linear 
programming  relaxation  of  random  LMAP  (5)  under  the  assumption  that  hyperedge  costs  (pir-id  are  iid 
uniform  on  [0,  1]: 

E[ Z5lf„l  < 

More  recently,  random  instances  of  the  Linear  MAP  (5)  have  been  analyzed  using  the  methods  of  sta¬ 
tistical  physics  by  Martin,  Mezard,  and  Rivoire  (2005)  in  an  attempt  to  generalize  the  corresponding 
conjecture  (12)  of  Mezard  and  Parisi  (1985)  for  the  random  LAP;  unfortunately,  the  approach  employed 


5 


by  the  authors  did  not  yield  an  explicit  expression  for  the  expected  optimal  cost  of  random  Linear  MAP 
(5). 

Krokhmal  et  al.  (2007)  have  established  a  limiting  value  of  the  expected  optimal  cost  E[Z^|„]  of  random 
Lineal-  MAP  (5)  with  absolutely  continuous  cost  distributions  F  that  admit  a  power  series  asymptotical 
expansion  in  the  neighborhood  of  the  left-end  point  F-1(0)  of  their  support  sets: 

lim  -E[ZdW]  =  F~\ 0)  =  inf  { f  |  F(t)  >0}.  (14) 

n  or  <i— *oo  n 

2  On  optimality  of  a  polynomial  algorithm  for  random  linear  multidi¬ 
mensional  assignment  problem 

The  Lineal'  Multidimensional  Assignment  Problem  (LMAP)  is  a  higher-dimensional  generalization  of 
the  well-known  two-dimensional,  or  Linear  Assignment  Problem  (LAP)  (see,  e.g.,  Papadimitrou  and 
Steiglitz,  1998;  Burkard  et  al.,  2009).  A  graph-theoretic  formulation  of  the  LAP  of  cardinality  n  presents 
it  as  finding  a  minimum-cost  perfect  matching  on  a  balanced  bipartite  graph  with  2 n  vertices,  provided 
that  the  cost  of  matching  is  defined  as  the  sum  of  edge  costs.  Similarly,  a  d -dimensional  LMAP  of 
cardinality  n  can  be  formulated  as  finding  a  perfect  matching  on  a  balanced  d -partite  d -uniform  hyper¬ 
graph  with  dn  vertices,  such  that  the  sum  of  the  costs  of  hyperedges  in  the  matching  is  minimized  (see, 
among  others,  Berge,  1989;  Bollobas,  1998  for  general  references  on  hypergraphs).  If  the  cost  of  hy¬ 
peredge  (i\ . id),  where  i\ . id  e  {1, . . . is  given  by  <f>iv~id ,  the  mathematical  programming 

formulation  of  the  LMAP  of  dimensionality  d  and  cardinality  n  reads  as 


i\  =  1 

ik  =  l,...,n,  (15) 

k  =  2 _ _  d  —  1, 

id  =  1 . n, 

ik  =  1 _ _  n,  k  =  1, _ d, 

where  Xir..,d  =  1  if  the  hyperedge  (i  i . id)  is  included  in  the  matching,  and  x, , ...  ld  =  0  otherwise. 

It  is  easy  to  see  that  a  feasible  solution  of  (15)  is  given  by  n  hyperedges  [i\  , ...  ,i\  ),  r  =  1 
such  that  in  each  dimension  k  the  set  \i^ . ij," *  [  is  a  permutation  of  {1 . n  j  .  Hence,  problem 

(15)  admits  the  following  geometric  interpretation:  given  a  d -dimensional  matrix  <t>  =  { r/j, , ... irl }  e  , 
find  such  a  permutation  of  its  “rows”  and  “columns”  in  all  dimensions  that  the  sum  of  the  diagonal 
elements  is  minimized,  and  thus  is  also  known  in  the  literature  as  “axial”  multidimensional,  or  multi¬ 
index  assignment  problem. 

The  LMAP  (15)  was  first  introduced  by  Pierskalla  (1968),  and  has  found  numerous  applications  in  the 
areas  of  data  association,  image  recognition,  multi-sensor  multi-target  tracking,  peer-to-peer  satellite 


n  n 

min  "■  J2  4>h~idxh~id 

*1  =  1  id  =  1 
n  n 

s-1'  £  =  !’ 

*2=1  id  —  1 

n  n  n  n 

E  E  •••  ^  xiv..id  =  l, 
i  1  =  1  ik- 1  =  1  *£+i  =  l  id  —  1 

n  n 

E  •••  L  xh-id  = 

6=1  id- 1  =  1 
xi\-id  e  {O’  (}’ 
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refueling,  and  so  on  (for  a  detailed  discussion  of  the  properties  and  applications  of  the  LMAP,  see,  for 
example,  Burkard  and  Cel  a,  1999b;  Burkard,  2002;  Burkard  et  ah,  2009). 

In  contrast  to  the  LAP  that  is  polynomially  solvable,  the  LMAP  with  d  >  3  is  generally  NP-complete, 
a  fact  first  established  by  Karp  (1972)  for  the  3-dimensional  (d  =  3)  assignment  problem  (15).  For 
a  discussion  of  approximation  properties  of  the  LMAP  (15),  see,  e.g.,  Spieksma  (2000);  in  particular, 
Crania  and  Spieksma  (1992)  have  shown  that  even  in  the  case  when  costs  0;ii-2?-3  of  the  3-dimensional 
LMAP  arc  decomposable,  i.e.,  </>;ii-2z3  =  r/,i;-2  +  <r/7_,l3  +  r/,i;-3,  there  is  no  polynomial  algorithm  that 
yields  an  £-approximate  solution  for  any  s  >  0. 

Exact  and  heuristic  algorithms  for  three-  and  higher-dimensional  LMAPs  were  proposed  in  Balas  and 
Saltzman  (1991);  Poore  (1994a,b);  Poore  and  Robertson  (1997);  Murphey  et  al.  (1998a,b);  Andrijich  and 
Caccetta  (2001),  and  others.  In  particular,  Balas  and  Saltzman  (1991)  introduced  a  number  of  heuris¬ 
tics  for  the  3-dimensional  LMAP;  Andrijich  and  Caccetta  (2001)  report  that  on  randomly  generated 
problems,  some  of  these  heuristics  yield  solutions  very  close  to  optimality.  These  observations  find  a 
theoretical  substantiation  in  Kravtsov  (2005),  who  has  demonstrated  that  if  the  assignment  costs  in  (15) 
are  iid  random  variables  from  a  discrete  distribution  satisfying  certain  properties,  a  simple  greedy  al¬ 
gorithm  produces  asymptotically  optimal  solutions  with  high  probability,  when  the  cardinality  n  of  the 
LMAP  (15)  increases  infinitely. 

In  this  work,  we  strengthen  and  generalize  the  results  of  Kravtsov  (2005),  by  showing  that  a  greedy 
algorithm  produces  e-approximate  solutions  of  random  LMAP  almost  surely  (a.  s.),  or,  in  other  words, 
that  the  cost  of  the  greedy  solution  converges  strongly  to  the  optimal  cost.  Further,  we  extend  the  analysis 
to  random  LMAPs  whose  costs  arc  continuously  distributed,  including  distributions  with  unbounded 
support  sets. 

Results  concerning  asymptotic  optimality  of  heuristic  algorithms  on  randomly  generated  problems  arc 
well  known  in  the  context  of  other  hard  combinatorial  optimization  problems,  such  as  the  Quadratic 
Assignment  Problem  (QAP),  which  is  also  known  to  be  NP-complete  and  non-approximable  (see,  among 
others,  Burkard  et  ah,  2009;  Krokhmal  and  Pardalos,  2009,  and  references  therein).  In  the  case  of  random 
QAP,  asymptotic  optimality  of  heuristic  solution  methods  is  a  manifestation  of  the  fact  that  for  instances 
of  random  QAP  large  enough,  all  its  feasible  solutions  are  asymptotically  optimal  (Burkard  and  Fincke, 
1982a).  Moreover,  an  entire  class  of  combinatorial  optimization  problems  exists  that  shares  this  property 
with  the  QAP  (Burkard  and  Fincke,  1985;  Szpankowski,  1995).  The  LMAP  (15),  however,  does  not 
belong  to  this  class;  recent  investigations  of  asymptotic  behavior  of  random  LMAPs  (Krokhmal  et  al., 
2007;  Krokhmal  and  Pardalos,  2011)  entail  that  only  a  vanishingly  small  fraction  the  feasible  set  of  a 
random  LMAP  is  e-optimal. 

2.1  A  greedy  algorithm  for  LMAP  with  discrete  iid  random  costs 

Algorithm  1  describes  the  greedy  heuristic  for  the  LMAP  that  is  in  the  focus  of  this  work.  The  heuristic 
starts  by  finding  the  smallest  hyperedge  cost  of  the  LMAP  (15),  which  we  denote  by  cf>. <n  (n,  and  then 

*1  "ld 

removing  from  the  cost  matrix  <L  the  cost  elements  0.  .(d  .  for  each  i  G  {1, . . . ,  d},  i.e.,  the  costs  of 

*1  ”’li  "'ld 

the  hyperedges  that  arc  not  feasible  with  respect  to  the  smallest-cost  hyperedge  (/['  * . z'^).  Then,  the 

procedure  is  repeated,  and  upon  finding  the  smallest  hyperedge  cost  cp.<2)  .<?)  in  the  reduced  cost  array 

h 

T>,  the  costs  of  hyperedges  that  are  infeasible  with  respect  to  the  hyperedges  (i[k^ . z^),  k  =  1,2, 

are  discarded,  and  so  on.  After  n  steps,  n  costs  r/> . a- >  .i/o  are  obtained,  which  have  the  property  that  for 
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each  l  —  1, . . . ,  d  the  indices  . J  are  all  different,  i.e.,  a  feasible  solution  of  (15)  is  found. 

Algorithm  1  A  greedy  heuristic  for  LMAP  (15) 
l:  input:  Cost  matrix  O  =  {<pil..-id  |  (/ 1, . . . ,  id)  €  {1, . . . ,  n}d}  e 
2:  initialize:  Z„  :=  0;  for  each  l  €  {1, . . . ,  d}  define  set  Mi  {1 

3:  for  k  :=  1  to  n  do 

4:  define  a  submatrix  e  of  the  cost  matrix  <f>  as 

<$>{k)  :=  {pii-ia  I  O't . id)  e  M  x  •  •  •  x  A/^} 

5:  find  the  smallest  element  0.ot)  .«)  of  the  submatrix 

0'P\  •  •  • ,  e  arg min  €  <I>(fc)} 

6:  let  Zn  \=  Zn  +  (j>  (k)  Ak) 

l\  '"v 

7:  for  each  €  e  {1, . . . ,  d}  update  the  set  Mi  :=  Mi\{i^  } 

8:  end  for 

9:  for  each  A:  e  {1 . «}  define  xaw  ■<*>  :=  1  andx,j.../rf  :=  0  for  all  other  (/ 1 _ ,id) 

l\  "'ld 

10:  output:  A  feasible  solution  x;-|  -  id  of  LMAP  (15)  and  its  cost  Zn 


Obviously,  the  described  greedy  heuristic  for  LMAP  runs  in  0(nd+l)  time.  The  next  lemma  provides  a 
foundation  for  the  subsequent  probabilistic  analysis  of  the  greedy  heuristic  and  is  a  strengthening  of  the 
corresponding  result  in  Kravtsov  (2005). 


Lemma  1.  Consider  a  set  Sn  of  cardinality  |5n|  =  \cn  whose  elements  are  iid  random  variables  dis¬ 
tributed  uniformly  over  pn  values  an  <  ...  <  bn.  Assume  that  pn  and  Kn  increase  with  n  such  that  the 
following  series  converges: 


1 2p"e 

n 


Kn 

Pn 


<  oo. 


Then,  for  n  sufficiently  large,  the  set  Sn  contains  the  minimum  element  an  almost  surely  (a.s.) 


Proof  To  verify  the  statement  of  the  lemma,  it  is  convenient  to  think  about  the  set  Sn  in  terms  of 
randomly  distributing  |5„|  =  Kn  different  objects  into  pn  boxes.  Then,  define  An  as  the  event  that  S„ 
contains  the  smallest  element  an ,  whence 

P {An}  <  P{at  least  one  box  is  empty}  =  1  —  P {B„},  (16) 

where  Bn  is  the  event  that  there  arc  no  empty  boxes,  for  which  it  holds  (see,  e.g..  Feller,  1968): 

_ Kn 

If  pn  and  Kn  increase  with  n  such  that  the  quantity  A„  =  pne  p><  is  bounded,  that  it  can  be  shown  that 
each  summand  in  the  above  sum  is  asymptotically  equal  to  {—Xn)l/i\,  whereby 


Thus,  for  n  sufficiently  large,  the  probability  that  the  set  Sn  does  not  contain  the  smallest  element  an  of 
the  distribution  can  be  bounded  as 

p {An}  <  1  -  p {Bn}  «  1  -  e~x"  =  Xn  +  0( X2n). 

Since  by  the  conditions  of  the  Lemma,  Xn  <  oo,  from  the  Borel-Cantelli  lemma  we  immediately 
have  that  P {An  i.o.}  =  0  P{An  ev.}  =1.  ■ 

Assuming  that  the  assignment  costs  of  the  LMAP  (15)  arc  positive,  a  feasible  solution  with  cost  Zn  is  an 
^-approximation  of  the  optimal  cost  Z*  of  LMAP  of  cardinality  n  if  it  satisfies 

Zn  5  Z*(l  +  e).  (17) 

The  next  theorem  establishes  the  conditions  on  the  discrete  distribution  of  assignment  costs  in  (15)  under 
which  the  greedy  algorithm  delivers  an  s- approximation  of  the  optimal  cost  of  the  LMAP,  or,  more 
precisely,  the  ratio  of  the  greedy  solution  cost  to  the  optimal  cost  approaches  unity  almost  surely. 

Theorem  1.  Consider  LMAP  (15)  with  d  >  3,  n  >  2,  whose  cost  coefficients  are  iid  random  variables 
distributed  uniformly  over  na  values 2  an  <  ...  <  bn,  where  an  >  0  and  a  >  0.  Then,  there  exists 
a  constant  M  >  0  such  that  the  greedy  algorithm  produces  a  solution  with  the  cost  Zn,  which  for 
sufficiently  large  n  satisfies 

—  -  l  <  M  ( —  -  \nl/d  n  a.  s.,  (18) 

Zn  \<*n  ) 

where  Z*  is  the  optimal  cost  of  the  LMAP. 

Proof.  From  the  description  of  the  greedy  heuristic,  it  follows  that  the  cost  of  the  feasible  solution  can 
be  represented  as 


n  n 

zn  =  y "i  #,■(*>...,•(*>  =  y (i9) 

k= 1  k= 1 

where  each  f]c  is  equal  to  the  smallest  element  of  a  submatrix  T(/l  ) : 

(f>k  =  min  {a  |  a  e  k  —  1 . n. 

In  general,  the  summation  in  (19)  contains  terms  (pk  that  arc  either  equal  to  the  smallest  element  an  of  the 
distribution,  or  exceed  it.  Let  Kn  denote  the  (random)  number  of  the  summands  in  (19)  that  are  greater 
than  an : 


Kn  =  \{k  I  4>k  >  an}\.  (20) 

Then,  noting  that  0  <  Kn  <  n,  the  optimal  cost  Z*  of  the  LMAP  (15)  and  the  cost  Zn  returned  by  the 
greedy  heuristic  can  be  bounded  as 

na„  <  Z*  <Zn  <  (n  -  K„)an  +  K„b„  =  na„  ( 1  +  —— — —  j  ,  (21) 

V  n  an  ) 

2Here  and  in  what  follows  we  omit  rounding  to  avoid  unnecessary  ramifications  in  exposition. 
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from  which  it  is  easy  to  see  that  Zn  is  an  ^-approximate  solution  of  (15)  by  means  of  the  approximation 
inequality  (17)  as  soon  as  Zn  <  na„  ( 1  +  e).  Thus,  for  some  en>  0  consider 

P{Z„  >  na„(l  +  sn ) }  <  P {Kn  >  ynn}, 

where  the  inequality  follows  from  (21)  provided  that  y„  is  chosen  as 

£n 

Yn  ~  7  7  7  ■ 

Dn/ttn  1 

Observe  that  if  Kn  >  y„n  holds,  then  there  exists  an  integer  v  €  {1 . n}  such  that  v  >  ynn  and 

the  corresponding  submatrix  of  $  with  vd  elements  does  not  contain  elements  equal  to  an.  Then,  from 
Lemma  4  it  follows  that  for  sufficiently  large  values  of  n , 

P  {Kn  >  ynn}  <  P{  set  of  size  ( ynn)d  does  not  contain  an  }  <  na  exp 

Choosing  the  parameter  e„  in  the  form 

=  (2 +  «)■/■<  (^t- 
\an 

we  have  that  for  values  of  n  large  enough, 

P  {Kn  >  ynn}  <  rr2, 

whence  expression  (18)  follows  by  the  Borel-Cantelli  lemma.  ■ 

Corollary  1.1.  If  a  <  d  in  (18)  and  the  ratio  hnlan  satisfies 

b?-  =  o(nl-a'd\n-l'dn),  n  »  1, 
an  s  ' 

then  for  sufficiently  large  n,  the  greedy  cost  Zn  is  an  e-approximation  of  the  optimal  cost  Z*  of  random 
LMAP  due  to  (17)  for  any  e  >  0  almost  surely.  Put  differently,  the  cost  ratio  Z„  /  Z  *  between  greedy  and 
optimal  solutions  converges  to  unity  a.  s.,  with  the  convergence  rate  given  by  (18). 

Remark  1.1.  The  intuition  behind  Lemma  4  and  Theorem  5  is  that  if  the  elements  of  the  cost  matrix  4> 
are  drawn  at  random  from  a  sufficiently  small  set  of  values,  then  at  each  step  of  the  greedy  heuristic,  the 
submatrix  will  contain  the  smallest  element  from  that  set  with  sufficiently  high  probability.  This 
observation  can  be  pressed  into  service  to  address  the  case  when  the  elements  of  the  cost  matrix  <t>  of  the 
LMAP  (15)  arc  continuous  iid  variables,  as  it  is  shown  next. 

2.2  Greedy  heuristic  for  LMAP  with  continuous  iid  costs 

In  (Krokhmal  et  ah,  2007;  Krokhmal  and  Pardalos,  201 1)  it  has  been  demonstrated  that  if  the  assignment 
costs  fii ! ...  i(]  in  LMAP  (15)  arc  iid  random  variables  with  a  continuous  distribution  F,  then  asymptotic 
behavior  of  the  optimal  value  of  random  LMAP  is  controlled  by  the  properties  of  the  distribution  F  in 
the  vicinity  of  the  left-end  point  F~x  (0)  of  the  support  set  of  the  distribution,  where 

F~]  (0)  =  inf  {t  |  F(t)  >  0}. 


1  I  n 


a/d—l 


ln1^  n , 
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In  view  of  that,  we  restrict  our  discussion  to  continuous  distributions  F  whose  support  sets  arc  bounded 
from  above, 

F-1(l)  <  +oo,  where  F_1(l)  =  sup{t  |  F(t)  <  1}. 

The  next  theorem  generalizes  the  results  of  the  previous  section  to  continuous  distributions. 

Theorem  2.  Consider  LMAP  (15)  with  d  >  3,  n  >  2,  whose  cost  coefficients  are  iid  random  variables 
with  a  continuous  distribution  F  that  has  a  bounded  support  [ a,b ],  where  a  >  0.  Then,  for  any  a  >  0, 
there  exists  a  constant  M  >  0  such  that  the  greedy  algorithm  produces  a  solution  with  cost  Zn,  which 
for  sufficiently  large  n  satisfies 

—  -  1  <  — — — — -  -  1  +  Mnald~x  ln1/rf  n  a.  s.,  (22) 

Z*  ~  a 

where  Z*  is  the  optimal  cost  of  the  LMAP. 

Proof.  For  a  continuous  distribution  F  on  [a,  b]  C  M,  define  the  sequence  {8n(k}},  k  —  0 . pn,  as 

4(0)  =  0,  8n(k)  -  -a  +  F~x  +  F(a  +  8n{k  -  1))^  ,  k  =  \,...,pn. 

The  intervals  2*  =  (a  +  8n(k  —  l),a  +  8n(k)],k  —  1, pn,  partition  the  set  (a,  b]  into  pn  equiprobable 
“bins”,  such  that  for  any  F -distributed  random  variable  X 

P {a  +  8n(k  —  1)  <  X  <  a  +  8n(k)}  =  — ,  k  —  1 . pn. 

Pn 

Then,  the  elements  of  the  cost  matrix  can  be  labeled  with  pn  different  labels,  in  accordance  to  the  “bin” 
Zfc  that  the  corresponding  cost  element  falls  into.  Obviously,  the  labels  are  independently  and  identically 
uniformly  distributed.  Therefore,  taking  into  account  that  the  elements  of  the  cost  matrix  T>  that  fall  into 
bin  X\  are  less  than  or  equal  to  a  +  8n  (1),  the  cost  Zn  of  the  greedy  solution  of  the  MAP  can  be  bounded 
as 


na  <  Z*  <  Zn  <  (n  -  Kn)(a  +  5„(1))  +  Knb.  (23) 

where  K„  is  equal  to  the  number  of  summands  in  the  cost  Zn  of  the  greedy  solution  that  do  not  fall  into 
the  bin  L\  ~  (a,  a  +  8n(  1)].  Then,  for  any  fixed  en  >  0  it  holds  that 

P{Z„(€>)  >  na(  1  +  £„)}  <  P { (n  -  K„)(a  +  4(1))  +  Knb  >  na{  1  +  £„)}  =  P  {Kn  >  nyn} , 

where 

_  s„  -  8n(\)/a 
b/a  —  1  —  8n(l)/a 

Similarly  to  the  arguments  of  Theorem  5,  Kn  >  nyn  holds  provided  that  there  exists  v  €  {1 . n}  such 

that  v  >  ny„  and  the  corresponding  submatrix  of  size  vd  does  not  contain  elements  from  the  interval  T\ , 
whence 

D  ,  „  ^  ^  (  (Ynn)d  l 

P  {Kn  >  ynn}  <  pn  exp  f - , 

I  Pn 
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Let  pn  —  na  for  some  a  >  0,  then,  choosing  en  as 


£11  — 


F~l{n~a) 


1  +  1- 
a 


b  F~'(n-a) 


a  \a  a 

—  zr-ir„-a\ 


(a  +  2)1/dna/d-1ln1/d , 


and  taking  into  account  that  8„(l)  =  F  1(n  “)  —  a  =  o(l),  n  1,  we  obtain  that 

P {Kn  >  ynn)  <  p~ 


~Va  =  n-2 


which  verifies  statement  (22)  of  the  Theorem  by  virtue  of  the  Borel-Cantelli  lemma.  ■ 

Corollary  2.1.  If  a  <  d,  it  follows  from  (22)  that  Zn  represents  an  e- optimal  solution  of  random  LMAP 
in  the  sense  (17)  for  any  e  >  0,  and  the  cost  ratio  Zn/Z*  converges  to  unity  a.  s.  It  is  natural  that  the 
value  of  the  parameter  a  6  (0,  d )  in  (22)  is  selected  based  on  the  properties  of  F~l  at  the  origin  so  as 
to  increase  the  rate  of  convergence.  In  particular,  if  in  some  neighborhood  of  0  the  inverse  F~l  of  the 
distribution  F  satisfies  for  some  v  >  0 

F~\u)  <  a  +  Luv,  L  >  0,  u  0+, 


then  there  exists  a  constant  M\  >  0  such  that 


—  —  1  <  Min('1+vd^~1~l  ln1^  n  a.s. 
Z* 


Next,  we  consider  the  case  of  a  continuous  distribution  F  with  support  of  the  form  (— oo,  b],  where  the 
following  bounds  on  the  optimal  cost  of  random  LMAP  play  a  key  role.  Namely,  as  shown  in  Krokhmal 
and  Pardalos  (2011),  the  optimal  value  Z*  of  random  LMAP  with  iid  cost  coefficients  whose  distribution 
has  a  support  unbounded  from  below,  satisfies  for  sufficiently  large  n 


nF 


~d- 1 


<  Z*  < nF~ 


3  Inn 

nd  1 


a.  s. 


(24) 


Expression  (24)  entails  that  when  support  of  F  is  unbounded  from  below,  F  1  (0)  =  — oo,  one  has  that 
Z*  <  0  a.  s.  for  large  enough  n.  Note  that  in  this  case  the  approximation  condition  (17)  takes  the  form 

Zn  <  Z*(l  -£),  £  >  0.  (25) 


Taking  into  account  (25),  the  following  statement  holds  regarding  the  quality  of  the  greedy  solution  to  a 
random  LMAP  (15). 


Theorem  3.  Consider  LMAP  (15)  with  d  >  3,  n  >  2  whose  cost  coefficients  are  iid  random  variables 
with  continuous  distribution  F  such  that  F-1(0)  =  — oo,  F-1(l)  <  oo.  Then,  the  greedy  algorithm 
produces  a  solution  with  cost  Zn  that  for  sufficiently  large  n  satisfies 


1-  —  < 

Z*  ~ 


d  —  1 


1  /d 


7-1 


+ 


,d- 1 


3  In  n 

nd  1 


-  1 


a.  s. 


where  Z*  is  the  optimal  cost  of  the  LMAP. 


(26) 
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Proof.  Similarly  to  the  proof  of  Theorem  2,  let  us  partition  the  semi-infinite  support  (—00,  b]  of  the 
distribution  F  into  pn  “bins”  {an(k  —  1),  c/n(k)\  such  that 

P{an(k  —  1)  <  X  <  a„(k)}  -  — ,  k  -  1  ,...,pn, 

Pn 

where  X  is  F-distributed  random  variable,  and  an(k )  is  defined  as 

MO)  =  -00,  an(k)  =  F_1  ^F(an(k  -  1))  +  ,  k  =  l,...,pn.  (27) 


Then,  similar  arguments  allow  us  to  construct  an  upper  bound  Zn  on  the  greedy  cost  Zn  in  the  form 

Zn<Zn  =  (n  -  Kn)F-\p-1 )  +  K„F~\  1), 

where  /V/c  is  the  number  of  summands  in  the  greedy  cost  Zn  that  do  not  fall  into  the  first  “bin” 
(—00,  F_1(p“1)]. 

Next,  observe  that  if  the  optimal  cost  Z*  of  the  LMAP  can  be  bounded  from  below  and  above,  e.g., 

7  <  7*  <  7 

then,  given  a  fixed  e  >  0,  the  greedy  solution  cost  Zn  satisfies  the  approximation  inequality  (25)  as  soon 
as  the  upper  bound  Zn  satisfies 


Zn  —  z„  <  —eZn 


In  view  of  the  bounds  (24)  on  the  optimal  cost  of  random  LMAP  due  to  (Krokhmal  and  Pardalos,  201 1) 
that  hold  almost  surely  for  large  enough  n ,  define 


Z„  =  n  F~ 


1 


3pn  In  n 

and  for  any  fixed  en  >  0  consider  the  probability 


,  Zn  =  nF  1  I  —  I  ,  where  p„  = 


Pn 


*d- 1 


3  In  n 


P  | Zn (T>)  —  Z„  >  SnZn |  —  P {Kn  >  nyn}, 


where 


Yn  =  \  F 


r-l 


1 


3  pn  In  n 


-  F~l  (  — )(l  +  £„))  ( F-J(  1)  -  F~l  1  1 


Pn 


Pn 


-1 


(28) 


Following  the  arguments  of  Theorems  5  and  2,  it  can  be  shown  that  for  sufficiently  large  values  of  n  the 
probability  in  (28)  satisfies 


P {Kn  >  ny„}  <  pn  exp 


(nyn)d\ 

Pn 


(29) 
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If  one  selects  the  parameter  en  in  the  form 


--1 


£n  — 


1 


3  pn  In  n 


F~l  (l On1 ) 


d-\\x/d  (  In (3 Inn)  \  x,d 


1  - 


(<r/  —  1)  In  n 


n 
F~l 


+ 


1 


3  Pn  In  77 


F~l  (P/71) 


- 1, 


inequality  (29)  implies  that  the  probability  in  (28)  is  bounded  as 

P{K„  >  ny„}  <  p~2  <  n~2 


for  all  large  enough  values  of  n ,  thereby  verifying  the  estimate  (26)  of  approximation  quality  of  the 
greedy  solution  by  means  of  the  Borel-Cantelli  lemma.  ■ 


Remark  3.1.  In  the  case  when  the  distribution  F  is  such  that 


F~l 


3  In  n  \ 

nd~x  J 


1, 


oo, 


(30) 


Theorem  3  asserts  that  the  solution  cost  produced  by  the  greedy  heuristic  is  an  £-approximation  of  the 
optimal  cost  of  random  LMAP,  for  any  s  >  0.  Condition  (30)  holds,  for  instance,  for  distributions 
F  whose  inverse  F~l  has  a  logarithmic  singularity  at  the  origin,  i.e.,  when  the  following  asymptotic 
representation  holds  in  the  vicinity  of  0: 


F~\u) 


-Co 


In*!, 


0+ 


for  some  cq  >  0,  /3  >  0. 


Continuous  distributions  that  satisfy  this  condition  and  whose  support  is  bounded  from  above  include, 
for  instance,  exponential  distribution  on  (— oo,  0],  truncated  normal  distribution  on  (— oo,  b]\ 


F(t)  —  e?l(_oo;o](0  +  l(o,oo)(0-  F(t) 


m 

<&(b) 


1(— OO,  b](t)  +  1(2, 

,oo  )(0, 


where  T(/ )  is  the  standard  normal  distribution  function.  Observe  that  in  the  case  when  the  inverse  F  1 
of  the  cost  distribution  F  has,  for  example,  a  power  singularity  at  the  origin,  i.e., 

F~l{u)  ~  —  cotr-*,  w  — ►  0+,  co,/6>0. 


the  ratio  in  (30)  is  unbounded  in  n,  hence  no  statement  can  be  inferred  from  Theorem  3  regarding  the 
quality  of  the  greedy  solution  in  this  case. 

Remark  3.2.  According  to  Krokhmal  and  Pardalos  (201 1),  the  asymptotic  behavior  of  the  optimal  cost 
of  random  LMAP  (15)  is  determined  completely  by  the  properties  of  the  distribution  function  F  in  the 
vicinity  of  the  left-end  point  of  its  support.  Nevertheless,  the  requirement  of  boundedness  from  above  of 
the  distribution’s  support,  F-1(l)  <  oo,  which  is  imposed  in  Theorems  5-3,  is  essential  for  estimating 
the  quality  of  the  greedy  solution,  as  it  allows  one  to  obtain  an  upper  bound  on  the  cost  of  the  solution 
produced  by  the  greedy  algorithm. 
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3  On  Hamming  distance  in  hypergraph  matching  problems 


The  concept  of  distance  between  feasible  solutions  of  optimization  problems  plays  an  important  role  in 
combinatorial  optimization.  For  instance,  it  is  widely  acknowledged  that  analysis  of  the  'problem's  fitness 
landscape  (Weinberger,  1990;  Stadler,  1996),  which  comprises  the  set  of  feasible  solutions,  their  fitness 
values  (costs)  / (•),  and  a  measure  d(- ,  •)  of  distance  between  solutions,  can  yield  useful  insights  into  the 
performance  and  tuning  of  exact  and  heuristic  algorithms. 

This  section  discusses  some  properties  of  the  Hamming  distance,  a  popular  distance  measure  d(- ,  ■)  in 
combinatorial  optimization,  in  application  to  problems  where  the  underlying  combinatorial  structure  is 
based  on  hypergraph  matchings,  which  generalize  the  well-known  class  of  combinatorial  problems  on 
bipartite  graph  matchings,  such  as  the  linear  assignment  problem  (LAP),  quadratic  assignment  problem 
(QAP),  etc.  (comprehensive  reviews  of  assignment  problems  can  be  found  in,  e.g.,  Pardalos  and  Pitsoulis, 
2000;  Burkard,  2002;  Burk  aid  et  al.,  2009). 


The  Hamming  distance  was  first  introduced  by  Hamming  (1950)  as  a  measure  of  errors  (or  substitutions) 
that  transform  one  string  of  a  binary  code  into  another,  and  since  then  has  found  applications  in  the  areas 
of  coding  theory,  information  theory,  cryptography,  combinatorial  optimization,  and  others  (Matsumoto 
et  ah,  2006).  Given  two  strings  of  equal  length  with  characters  from  any  alphabet  (not  necessarily  binary), 
the  Hamming  distance  between  them  is  usually  defined  as  the  number  of  positions  in  which  these  strings 
disagree.  The  landscape  structure  of  many  combinatorial  optimization  problems  can  be  investigated 
using  the  Hamming  distance.  In  particular,  the  Hamming  distance  defined  on  the  set  of  permutations  of 
a  given  length  was  applied  to  study  the  fitness  landscape  of  the  quadratic  assignment  problem  (Merz  and 
Freisleben,  2000).  The  feasible  set  of  the  QAP  constitutes  a  special  case  of  (1)  with  d  —  2,  and  by  virtue 
of  (10)  a  feasible  solution  of  the  QAP  can  be  represented  in  the  form  p  =  ((,  n),  whereby  the  Hamming 
distance  between  two  solutions  pn  =  (i ,  717 )  and  pj  =  (t ,  nj )  is  equal  to  the  number  of  positions  in 
which  nj  (k)  jtj(k),k  =  1 


It  is  easy  to  see  that  such  a  definition  does  not  apply  in  the  case  of  multidimensional  assignment  prob¬ 
lems  (3)  with  d  dimensions  and  n  elements  per  dimension,  where  a  feasible  solution  is  generally  an 
unordered  collection  of  n  strings  of  length  d.  Thus,  for  hypergraph  matching  problems  the  Ham¬ 


ming  distance  can  be  defined  in  terms  of  the  minimum  number  of  positions  in  which  two  matchings 


Mi 


=  {Of, 


.*5°). 


o-r. 


ijl))\  and  pj  =  lO'i 


;(D 


41}) . of, 


.j(dn))}  differ 


from  each  other.  For  example,  the  Hamming  distance  between  the  following  two  feasible  solutions  of  a 
d  =  4  and  n  —  3  multidimensional  assignment  problem: 


/  1 1 1 1  \  /  1222  \ 
Pi  =  2222  ,  p2  =  2111  , 

V  3333  /  V  3333  ) 


is  equal  to  2,  but  not  6.  In  general,  let  us  define  the  number  of  elements  by  which  the  k-th  hyperedge  in 
Pi  differs  from  the  l- th  hyperedge  in  pj  as 


A kl  =  (1 


:(k) 


f ) 


Of. 


■jf)  II  = 


r=  1 


(Ar)  •(£), 
r  Jr 


(31) 


where  Sp  is  the  negation  of  Kronecker’s  delta,  Sp  =  1  —  Sp.  One  evidently  has  =  A/;/c,  and 
0  <  A kt  <  d.  Then,  the  Hamming  distance  <:///(■ ,  •)  between  hypergraph  matchings  is  defined  as  the 
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optimal  value  of  the  following  LAP: 


n 

=  WlM  ~  dj  II  =  min  Y  Ak^{k),  (32) 

7Z€.\\n  , 

K  =  1 

where  Iln  is  the  set  of  all  permutations  n  :  [n]  i->  [n].  Below  we  discuss  some  properties  of  the 
introduced  Hamming  distance  measure  djjf ,  •)  as  an  optimal  value  of  the  linear  assignment  problem 
(32). 


3.1  Diameter  of  the  feasible  set 

Given  a  distance  measure  d{- ,  •)  in  a  combinatorial  optimization  problem,  the  diameter  of  the  feasible 
set  can  be  defined  as  the  maximum  distance  between  two  feasible  solutions.  The  Hamming  diameter  D 
of  the  multidimensional  assignment  problem  (3)  is  then  defined  as  the  maximum  value  of  the  Hamming 
distance  (32): 


D  =  max  dn(Ui.  A/)-  (33) 

Then,  we  have  the  following  simple  observation: 

Proposition  1.  The  Hamming  diameter  D  of  the  feasible  set  of  a  multidimensional  assignment  problem 
(3)  of  dimensionality  d  and  cardinality  n  satisfies 

D  <  n(d  —  1),  (34) 

with  the  equality  attained  for  problems  with  n  >  d. 

Proof  Inequality  (34)  for  D  follows  trivially  from  the  permutation  representation  (10)  of  a  feasible 
solution  of  the  MAP  (3).  To  show  that  this  bound  is  exact  for  problems  with  n  >  d,  consider  without 

loss  of  generality  the  distance  between  the  trivial  solution  p\  =  {(1, . . . ,  1),  (2 . 2) . ( n , . . . ,  n )} 

given  by  identity  permutations  m . Jid-\  —  1  in  (10)'  and  a  solution  jl  =  (ifi, . . . ,  if^),  where 

n i  =  l,  and  =  a(%_j),  with  cr  being  the  forward  cyclical  permutation.  Since  there  are  n  >  d 
different  cyclical  permutations  of  [[«]  (including  the  identity  permutation),  we  have  that  for  the  feasible 
solutions  pi  and  jl  the  quantities  defined  by  (31)  satisfty 

d 

Aki  =  X]  4,7 r,  (€)  =  d  -  1, 

r=  1 


whence  dn(pi^  fi)  —  n(d  —  1).  ■ 

3.2  Expected  distance  to  global  optimum 

Many  heuristic  solution  algorithms  for  combinatorial  problems  rely,  at  least  partly,  on  local  search.  In 
the  context  of  assignment  problems,  local  search  is  often  conducted  by  permuting  several  (usually,  two) 
positions  in  the  current  solution.  In  this  respect,  it  is  of  interest  to  estimate  the  number  of  permutations 
that  are  necessary  to  reach  a  global  optimum  from  a  current  feasible  solution,  or,  in  other  words,  the 
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Hamming  distance  from  the  given  solution  to  the  optimal  one.  To  this  end,  it  is  convenient  to  assume 
that  the  matching  cost  function  <b(/T)  in  (3)  has  the  form 

<W  =  LI  <t>eh eim  >  (35) 

^i\  vs  ^im 

where  JJ  is  an  operator  defined  over  some  set  of  cost  elements  {(/)}  indexed  by  the  hyperedges  of 
If  the  cost  elements  f  are  independent  identically  distributed  (iid)  random  variables  with  a  continuous 
distribution,  then  problem  (3)  has  a  unique  global  minimum  (since  the  costs  of  all  feasible  solutions 
arc  different  almost  surely),  and  the  location  of  the  optimal  matching  is  uniformly  distributed  over  the 
feasible  set. 

Let  H  denote  the  expected  Hamming  distance  between  the  (unique)  global  minimum  and  a  given  feasible 
solution  of  a  random  MAP;  then  the  expected  value  of  H  can  be  computed  as 

n\cl~ 1 

E [H]  =  E[E[//  |  Y]]=  EiH  I  F]  P(F  =  kh  (36) 

k=  1 


where  the  random  variable  Y  takes  value  k  ( I  <  k  <  n\d~l)  if  k- th  feasible  solution  is  the  global 
minimum,  and  zero  otherwise,  and  E [H  \  Y]  is  the  conditional  expectation  of  Hamming  distance  to  the 
global  optimum  given  its  location.  Since  0  <  H  <  n(d  —  l)  per  Proposition  1 ,  one  has 

n(d— 1)  n{d— 1) 

E  [H\Y]=  pP{H=p\Y}=  J2  P  -JcT-v'  (37) 

p= 0  p= 0  ”• 

where  Np  is  the  number  of  feasible  solutions  located  at  the  Hamming  distance  p  from  the  given  solution. 
If  the  assignment  costs  {0}  in  (35)  are  iid  continuous,  then 


E[tf] 


n(d  —  1) 

E  r 


P= 0 


Np 

n  ld~1 


(38) 


Again,  given  that  the  location  of  the  global  minimum  is  uniformly  distributed,  the  expression  for  E  [  //  J 
above  can  be  interpreted  as  the  expected  diameter  of  the  feasible  set  of  (3).  The  expected  Hamming 
distance  E  [  //  ]  for  any  distribution  of  assignment  costs  can  be  computed  explicitly  in  the  special  case 
n  =  2. 


Proposition  2.  Consider  a  hypergraph  matching  problem  (3)  with  cardinality  parameter  n  =  2,  and 
cost  function  (35)  where  the  cost  elements  {</>}  are  iid  random  variables  with  a  continuous  distribution. 
Then,  the  expected  Hamming  distance  E [H]from  a  given  solution  to  the  global  optimum  is  equal  to 


E  [H] 


Ie'KI 

1  2d  m°d  ^ 

d 

2 

(39) 


where  d  mod  2  is  the  remainder  on  division  of  d  by  2. 
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Proof.  Since  in  all  dimensions  of  an  n  —  2  MAP  there  arc  only  two  elements,  then  Hamming  distance 
between  two  feasible  solutions  is  always  even:  dnlpi,  ll  j )  =  0, 2,  4, ... ,  and  equals  to  the  number  of 
dimensions  in  which  the  two  solutions  differ.  It  is  easy  to  see  that  permuting  elements  in  all  d  dimensions 

of  a  solution  leaves  it  unchanged;  similarly,  making  d  —  1  permutations  in  dimensions  2,3 . d  is  the 

same  as  permuting  dimension  1,  and  so  on.  Hence,  the  Hamming  distance  between  two  solutions  of  an 
n  —  2  MAP  may  take  values 

H  =  2s,  5  =  0,1 . I  4  I  , 


where  s  is  the  number  of  dimensions  permuted.  The  number  of  feasible  solutions  that  are  at  a 
distance  H  =  2s  from  the  given  solution  is  given  by  (d)  for  s  =  1, . . . ,  \_d/2\  —  1.  A  special  case  arises 
when  d  is  even  and  s  —  [d / 2\  =  d/2,  where  the  number  of  solutions  that  arc  at  a  distance  H  =  2s  =  d 

~  ,d\A  ,  ■  (  1212  \  J  (  2121 

1  *  niiptn  Qvmmptrv  tnr  inQtonpp  I  I  unn  I 

1212 


from  a  given  solution  is  equal  to  1(1  due  to  symmetry;  for  instance,  , 

z  s  y  2121 

the  same  solution.  Observe  that  this  does  not  occur  when  s  =  \_d /2\  and  d  is  odd. 


represent 


In  the  general  case,  computation  of  the  number  Np  of  feasible  solutions  that  are  at  a  given  distance  p 
from  a  specified  solution  presents  significant  challenges  due  to  the  combinatorial  nature  of  its  definition 
via  a  solution  of  the  LAP  (32).  An  upper  bound  on  Np  can  be  obtained  by  ignoring  these  combinatorial 
considerations. 

Proposition  3.  In  a  hypergraph  matching  problem  (3)  with  d  <  n,  the  number  Np  of  feasible  solutions 
located  at  a  distance  p,  2  <  p  <  n(d  —  l),  from  a  given  solution,  is  bounded  from  above  as 

Np<  n  (n where  (^V!-  (40) 

On  \  /  *=t  \Pl)  7'=o  W 

and  the  summation  is  over  all  vectors  (p\, . . . ,  pp)  such  that  Yi  =  \  Pi  =  P,  0  <  pi  <  n,  pi  f  1. 


Proof  Expression  (40)  is  obtained  by  selecting  k  out  of  d  dimensions  in  representation  (9)  of  a  feasible 

solution  of  problem  (3),  and  permuting  /?,  out  of  n  elements  in  each  dimension,  where  the  ( p\ . pp) 

satisfy  the  above  conditions.  To  make  sure  that  representation  matix  of  the  permuted  solution  has  the 
largest  possible  number  of  entries  different  from  the  corresponding  elements  of  the  original  matrix,  we 
count  only  those  permutations  that  do  not  leave  any  of  the  pt  permuted  elements  unchanged.  The  number 
of  such  permutations,  or  derangements,  is  given  by  D(pi)  in  (40)  (see,  e.g.,  Stanley,  1986).  ■ 


Since  bound  (40)  ignores  the  combinatorial  structure  imposed  by  (32),  it  is  reasonable  to  expect  that 
it  can  be  rather  loose,  especially  for  larger  values  of  p,  when  the  combinatorial  effects  will  be  most 
prominent.  However,  for  smaller  values  of  p,  pertinent  to  solution  algorithms,  the  combinatorial  effects 
of  definition  (32)  should  be  lessened,  making  the  bound  (40)  tighter.  These  conclusions  are  supported  by 
the  numerical  studies  of  Hamming  distance  in  MAPs  that  arc  presented  next. 


3.3  Numerical  results 

In  this  section  we  report  the  results  of  computational  experiments  conducted  on  determining  the  expected 
Hamming  distance  to  the  global  optimum  (or,  equivalently,  to  the  given  feasible  solution)  E[i/],  and  the 
distribution  of  the  numbers  Np  of  feasible  solutions  located  at  a  Hamming  distance  p  from  the  given 


18 


solution.  The  computational  studies  involved  exhaustive  enumeration  of  the  feasible  set  of  a  multidimen¬ 
sional  assignment  problem  (3)  and  computing  the  Hamming  distance  (32)  between  each  pair  of  feasible 
solutions  using  the  shortest  augmenting  path  algorithm  for  dense  LAPs  (Jonker  and  Volgenant,  1987). 

Due  to  the  enumerative  nature  of  this  case  study,  only  small  size  MAPs  (d  =  3 . 8 ,  n  =  2, ...  ,5) 

were  examined. 

Figure  1  displays  the  expected  Hamming  distance  E[/7]  to  the  global  minimum  solution  of  the  MAP  (3). 
Interestingly,  E[/7]  exhibits  practically  linear  growth  when  one  of  the  parameters  d  or  n  is  fixed  and  the 
other  increases.  This  observation  supports  our  interpretation  of  E [ 77 J  as  the  expected  diameter  of  the 
feasible  set  of  (3). 
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Figure  1:  Expected  Hamming  distance  from  a  given  solution  (global  optimum)  in  multidimensional 
assignment  problem  (3)  for  fixed  dimensionality  cl  or  cardinality  n. 


Figure  2  displays  typical  distributions  of  the  number  Np  of  feasible  solutions  that  are  located  exactly  at  a 
Hamming  distance  p  from  the  given  solution.  In  particular,  the  distributions  of  Np  are  generally  skewed 
towards  larger  values  of  p.  Table  1  presents  a  comparison  of  the  numerically  determined  values  of  Np 
and  the  corresponding  upper  bound  as  given  by  Proposition  (3).  As  expected,  the  upper  bound  (40)  is 
very  loose  for  larger  values  of  the  Hamming  distance  p.  when  the  combinatorial  properties  of  definition 
(32)  arc  dominant,  and  is  quite  tight  for  smaller  values  of  p,  when  the  combinatorial  effects  of  (32)  arc 
less  pronounced. 
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Figure  2:  Number  Np  of  feasible  solutions  of  a  multidimensional  assignment  problem  (3)  located  at  a 
Hamming  distance  p  from  a  given  solution  (global  optimum). 
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Table  1:  Illustration  of  the  tightness  of  the  upper  bound  (40)  for  the  number  Np  of  feasible  solutions  of 
MAP  ( d  =  4,  n  =  4)  located  at  a  distance  p  from  the  given  solution. 


p 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

Bound 

24 

32 

252 

576 

1896 

4320 

6390 

12416 

12744 

7776 

729 

NP 

24 

32 

234 

576 

1488 

2736 

4293 

3552 

864 

0 
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4  High-quality  Solution  Sets  in  Randomized  Multidimensional  Assign¬ 
ment  Problems 

In  this  section  two  methods  will  be  described  that  can  be  used  to  obtain  mathematically  proven  high- 
quality  solutions  for  MAPs  with  large  cardinality,  or  large  dimensionality.  These  methods  utilize  the 
concept  of  index  graph  of  the  underlying  hypergraph  of  the  problem. 

4.1  Random  Linear  MAPs  of  Large  Cardinality 

In  the  case  when  the  cost  <t>  of  hypergraph  matching  is  a  linear  function  of  hyperedges’  costs,  i.e.,  for 
MAPs  with  lineal-  objectives,  a  useful  tool  for  constructing  high  quality  solutions  for  instances  with  large 
cardinality  (n  1)  is  the  so-called  index  graph.  The  index  graph  is  related  to  the  concept  of  line  graph, 
in  that  the  vertices  of  the  index  graph  represent  the  hyperedges  of  the  hypergraph. 

Namely,  by  indexing  each  vertex  of  the  index  graph  Q*  =  (V*.  £*)  by  (i\ . id )  e  {1 . n}  , 

identically  to  the  corresponding  hyperedge  of  TLd\n-  the  set  of  vertices  V*  can  be  partitioned  into  n 
subsets  V*,  also  called  levels ,  which  contain  vertices  whose  first  index  is  equal  to  k: 

n 

V*=  (J  •  V*  =  {(k,  i2,  ...,id)\h . id  6  {1,  •  •  • ,  n}}. 

k=  1 

For  any  two  vertices  i,  j  e  V*,  an  edge  (/,  j )  exists  in  Q* ,  (/,  j )  e  £*,  if  and  only  if  the  corresponding 
hyperedges  of  TLd\n  do  not  have  common  nodes.  In  other  words, 

£*  =  {O'.  7)  |  i  =  O't . id)J  =  O't . *'rf) :  ik  ±  Jk ,  k  =  1 . «}. 

Then,  that  it  is  easy  to  see  that  Q*  has  the  following  properties. 

Lemma  2.  Consider  a  complete,  d -partite,  n-uniform  hypergraph  Hd \n  —  (V,5),  where  \£\  =  nd ,  and 
V  =  U]fc=i  such  that  Vk  FI  V/  =  0,  k  ^  l  and  \V^\  —  n,  k  =  \ d.  Then,  the  index  graph 
Q*  =  (■ V*,£*)of'Hd\n  satisfies: 

1.  Q*  is  n-partite,  namely  V*  =  U)t=i  V *  fl  V *  =  &fori  ^  j ,  where  each  is  an  independent 

set  in  V*.'  for  any  i,  j  e  one  lias  (i.  j )  f  £* 

2.  |  |  =  nd~l 2 3  for  each  k  =  1 _ _  n 

3.  The  set  of  perfect  matchings  in  TLd\n  is  isomorphic  to  the  set  of  n -cliques  in  Q* ,  i.e.,  each  perfect 
matching  in  TLd\n  corresponds  uniquely  to  a  (maximum)  clique  of  size  n  in  Q*. 
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Figure  3:  A  3-uniform  3-partite  hypergraph  'Hj\n  and  its  index  graph  Q* .  The  vertices  of  Q*  shaded  in 
grey  represent  a  clique  (or,  equivalently,  a  perfect  matching  on  %d\n)- 


Let  us  denote  by  Q*  (a,,)  the  induced  subgraph  of  the  index  graph  Q*  obtained  by  randomly  selecting  an 
vertices  from  each  level  Vji  of  Q* ,  and  also  define  N (an )  to  be  the  number  of  cliques  in  Q*(an),  then 
based  on  the  following  lemma  Krokhmal  et  al.  (2007)  one  can  select  an  in  such  a  way  that  Q*(an)  is 
expected  to  contain  at  least  one  n -clique: 

Lemma  3.  The  subgraph  Q*(an)  is  expected  to  contain  at  least  one  n-clique,  or  a  perfect  matching  on 
'Hd\n  (i-e.,  E[N(an)]  >  1)  when  <yn  is  equal  to 


Ctn 


n 


d- 1 


/?! 


d- 1 
n 


(41) 


In  the  case  when  the  cost  coefficients  4>iv~id  of  MAP  with  linear  or  bottleneck  objective  are  drawn 
independently  from  a  given  probability  distribution,  Lemma  3  can  be  used  to  construct  high  quality 
solutions.  The  approach  is  to  create  the  subgraph  Q*mn(an),  also  called  the  a -set,  from  the  index  graph 
Q*  of  the  MAP  by  selecting  an  nodes  with  the  smallest  cost  coefficients  from  each  partition  (level)  of  Q* . 
If  the  costs  of  the  hyperedges  of  %d\n,  or,  equivalently,  vertices  of  Q* ,  are  identically  and  independently 
distributed,  the  a-sct  is  expected  to  contain  at  least  one  clique,  which  represents  a  perfect  matching  in 
the  hypergraph  TLd\n-  It  should  be  noted  that  since  the  a-sct  is  created  from  the  nodes  with  the  smallest 
cost  coefficients,  if  a  clique  exists  in  the  a-set,  the  resulting  cost  of  the  perfect  matching  is  expected  to 
be  close  to  the  optimal  solution  of  the  MAP 

Importantly,  when  the  cardinality  n  of  the  MAP  increases,  the  size  of  the  subgraph  Q*(an)  or  (7*in(a?n) 
grows  only  as  0(n),  as  evidenced  by  the  following  observation: 

Lemma  4 .  If  d  is  fixed  and  n  — oo,  then  an  monotonically  approaches  a  finite  limit: 


an  y  a  :=  \ed  1_|  as  n  /■  oo.  (42) 

Corollary  3.1.  In  the  case  of  randomized  MAP  of  large  enough  cardinality  n  3>  1  the  subset  (/*  m 
expected  to  contain  a  high-quality  feasible  solution  of  the  MAP  can  simply  be  chosen  as  QTfioe),  where 
a  is  given  by  (42). 

Observe  that  using  the  a-sct  g * in  (a)  for  construction  of  a  low-cost  feasible  solution  to  randomized  MAP 
with  linear  or  bottleneck  objectives  may  prove  to  be  a  challenging  task,  since  it  is  equivalent  to  finding  an 
n -clique  in  an  -partite  graph;  moreover,  the  graph  £/*jn(a)  is  only  expected  to  contain  a  single  n -clique 
(feasible  solution).  The  following  variation  of  Lemma  3  allows  for  constructing  a  subgraph  of  Q*  that 
contains  exponentially  many  feasible  solutions: 
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Lemma  5.  Consider  the  index  graph  Q*  of  the  underlying  hypergraph  jiq\n  °fa  randomized  MAP,  and 
let 


Pn 


2 


n 


d- 1 


n ! 


d- 1 
n 


(43) 


Then,  the  subgraph  G*(Pn)  is  expected  to  contain  2n  n-cliques,  or,  equivalently,  perfect  matching  on 

'Ud\n- 


Proof.  The  statement  of  the  lemma  is  easy  to  obtain  by  regarding  the  feasible  solutions  of  the  MAP 
as  paths  that  contain  exactly  one  vertex  in  each  of  the  n  “levels”  V*, . . . ,  V*  of  the  index  graph  G*  ■ 
Namely,  let  us  call  a  path  connecting  the  vertices  (1,  i^\  . . . ,  i^)  e  V*,  (2,  z®  > . . . ,  z  ®)  e  Vj, 

(n,  i ^ . i^)  €  V*  feasible  if  {z^,  ij^\ ....  i^}  is  a  permutation  of  the  set  {1 . n}  for  every 

k  —  2 . cl .  Note  that  from  the  definition  of  the  index  graph  G*  it  follows  that  a  path  is  feasible  if  and 

only  if  the  vertices  it  connects  form  an  zz-clique  in  G*  ■  Next,  observe  that  a  path  in  G*  chosen  at  random 

is  feasible  with  the  probability  (fpr  j  ,  since  one  can  construct  /jb(4-i)  different  (not  necessarily 
feasible)  paths  in  G*  ■  Then,  if  we  randomly  select  fn  vertices  from  each  set  V*  in  such  a  way  that  out  of 
the  (fn)n  paths  spanned  by  G*(Pn)  at  least  2"  are  feasible,  the  value  of  fn  must  satisfy: 

from  which  it  follows  immediately  that  fn  must  satisfy  (43).  ■ 

Corollary  3.2.  If  d  is  fixed  and  n  ->  oo,  then  fn  monotonically  approaches  a  finite  limit: 

Pn  /  f  :=  \2ed~1]  as  n  /■  oo. 

Remark  1.  Since  the  value  of  the  parameter  pn  (43)  is  close  to  the  double  of  the  parameter  an  (41), 
the  subgraph  G*mn(Pn),  constructed  from  selecting  /)„  nodes  with  the  smallest  cost  coefficients  from  each 
partition  (level)  ofG*  will  be  called  the  “2a -set”,  or  G* (2a ). 


Following  Krokhmal  and  Pardalos  (2011),  the  costs  of  feasible  solutions  of  randomized  MAPs  with 
lineal-  or  bottleneck  objectives  that  are  contained  in  the  a-  or  2a-sets  can  be  shown  to  satisfy: 

Lemma  6.  Consider  a  randomized  MAP  with  linear  or  bottleneck  objectives,  whose  cost  coefficients 
are  iid  random  variables  from  a  continuous  distribution  F  with  a  finite  left  endpoint  of  the  support, 
F~\  0)  >  — oo.  Then,  for  a  fixed  cl  >  3  and  large  enough  values  of  n,  if  the  subset  (7*in  (a)  (or, 
respectively,  G*ulfP))  contains  a  feasible  solution  of  the  MAP,  the  cost  Zn  of  this  solution  satisfies 

l"-')rl(«)+rl(^r)^"frlp)'  ”»>•  <44> 


in  the  case  of  MAP  with  linear  objective  (5),  while  in  the  case  of  MAP  with  bottleneck  objective  (7)  the 
cost  Wn  of  such  a  solution  satisfies 


t 


,d  —  1 


<Wn<F 


-t 


3  lnzz 

n  d  1 


n  »  1. 


(45) 
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4.2  Random  MAPs  of  Large  Dimensionality 


In  cases  where  the  cardinality  of  the  MAP  is  fixed,  and  its  dimensionality  is  large,  d  1,  the  approach 
described  in  section  4.1  based  on  the  construction  of  a-  or  2a-subsct  of  the  index  graph  Q*  of  the  MAP 
is  not  well  suited,  since  in  this  case  the  size  of  G*(a)  grows  exponentially  in  cl. 


However,  the  index  graph  Q*  of  the  underlying  hypergraph  'Hj\n  of  the  MAP  can  still  be  utilized  to 
construct  high-quality  solutions  of  large-dimensionality  randomized  MAPs. 


Let  us 


call  two 
;  H)\ 


matchings  /jl,  =  {(7, 


(i) 


j  on  the  hypergraph  'Hd\n  disjoint  if 


in) 


(n) 


i^)}  and  m 


An) 


( if\ 


Ak ) 


)  +  0, 


;(0 


;f) 


for  all  1  <  k.  I  <  n. 


or,  in  other  words,  if  m  and  ji  j  do  not  have  any  common  hyperedges.  It  is  easy  to  see  that  if  the 
cost  coefficients  of  randomized  MAPs  arc  iid  random  variables,  then  the  costs  of  the  feasible  solutions 
corresponding  to  the  disjoint  matchings  arc  also  independent  and  identically  distributed. 

Next,  we  show  how  the  index  graph  Q*  of  the  MAP  can  be  used  to  construct  exactly  nd~x  disjoint 
solutions  whose  costs  arc  iid  random  variables.  First,  recalling  the  interpretation  of  feasible  MAP  solu¬ 
tions  as  paths  in  the  index  graph  Q* ,  we  observe  that  disjoint  solutions  of  MAP,  or,  equivalently,  disjoint 
matchings  on  arc  represented  by  disjoint  paths  in  Q*  that  do  not  have  common  vertices. 

Note  that  since  each  level  Vj*  of  Q*  contains  exactly  nd~ 1  vertices  (see  Lemma  2),  there  may  be  no  set 
of  disjoint  paths  with  more  than  nd~l  elements. 


On  the  other  hand,  recall  that  a  (feasible)  path  Q*  can  be  described  as  a  set  of  n  vectors 

^  =  {(i1(1),...,i‘,)) . of . f», 

such  that  { i  ^ . i^}  is  a  permutation  of  the  set  {1 . n }  for  each  k  =  1 . d.  Then,  for  any 

given  vertex  v ( 1  ^  =  (1,  f1' . z  9  ^ )  €  V*,  let  us  construct  a  feasible  path  containing  v ( 1  ^  in  the  form 


!(l.4° . ^1)),(2.«f) . 


)}. 


where  for  k  —  2, 


,  d  and  r  =  2 . n 

f  ;('—!) 


= 


+  L 

1. 


if  i 


Ar- 1) 


=  1 . n  —  1 , 


•  j:  .  (r—  1) 

if  i,  =  n. 


(46) 


In  other  words,  {i  5 1  , . . . ,  ij"'  |  is  a  forward  cyclic  permutation  of  the  set  { 1 }  for  any  k  =  2, ...  ,d . 

Applying  (46)  to  each  of  the  nd -1  vertices  (1  ,i^ . z’9^)  e  V*,  we  obtain  nd~l  feasible  paths 

(matchings  on  jjd\n)  that  arc  mutually  disjoint,  since  (46)  defines  a  bijective  mapping  between  any  vertex 

(hyperedge)  (k.  i^ . i^'  )  from  the  set  V?,  k  —  2, . . .  ,n,  and  the  corresponding  vertex  (hyperedge) 

U(1)  €  V*. 


Then,  if  hyperedge  costs  4>i\-id  in  the  linear  or  bottleneck  MAPs  (5)  and  (7)  are  stochastically  inde¬ 
pendent,  the  costs  4>(/Xi), _ ^(p,nd-i)  of  the  z?^-1  disjoint  matchings  pt\, _ And~x  defined  by  (46) 

arc  also  independent,  as  they  do  not  contain  any  common  elements  <piv-id.  Given  that  the  optimal  solu¬ 
tion  cost  Z*d  n  (respectively,  Wd  )  of  randomized  linear  (respectively,  bottleneck)  MAP  does  not  exceed 
the  costs  . . ,  &(pLnd-\ )  of  the  disjoint  solutions  described  by  (46),  the  following  bound  on  the 

optimal  cost  of  linear  or  bottleneck  randomized  MAP  can  be  established: 


23 


Lemma  7.  The  optimal  costs  Z*d  ,  Wf  n  of  random  MAPs  with  linear  or  bottleneck  objectives  (5),  (7), 
where  cost  coefficients  are  iid  random  variables,  satisfy 


7*  <  yY 

Xd,n  - 


W*  <  Ymax 
Wd,n  -  XV.nd- 1’ 


(47) 


where  Xp,  Xfax  ( i  =  1, . . . ,  nd~l)  are  iid  random  variables  with  distributions  F^-,max  that  are  deter¬ 
mined  by  the  form  of  the  corresponding  objective  function,  and  X  j  :/c  denotes  the  minimum-order  statistic 
among  k  iid  random  variables. 


Remark  2.  Inequalities  in  (47)  are  tight:  namely,  in  the  special  case  of  random  MAPs  with  n  =  2,  all 
of  the  n \d~ 1  =  2d~ 1  feasible  solutions  are  stochastically  independent  Grundel  et  al.  (2007),  whereby 
equalities  hold  in  (47). 


As  shown  in  Krokhmal  and  Pardalos  (201 1),  the  following  quality  guarantee  on  the  minimum  cost  of  the 
nd~x  disjoint  solutions  (46)  of  linear  and  bottleneck  MAPs  can  be  established: 

Xpnd_x  <  nF-1  (//~^)  ,  <  F~l  (n~^)  ,  d  »  1, 

where  F~l  is  the  inverse  of  the  distribution  function  F  of  the  cost  coefficients  <f>iv..id.  This  observation 
allows  for  constructing  high-quality  solutions  of  randomized  linear  and  bottleneck  MAPs  by  searching 
the  set  of  disjoint  feasible  solutions  as  defined  by  (46). 


4.3  Numerical  Results 

Sections  4.1  and  4.2  introduced  two  methods  of  solving  randomized  instances  of  MAPs  by  constructing 
subsets  (neighborhoods)  of  the  feasible  set  of  the  problem  that  arc  guaranteed  to  contain  high-quality 
solutions  whose  costs  approach  optimality  when  the  problem  size  (n  -v  oo,  or,  respectively,  d  ->  oo) 
increases.  In  this  section  we  investigate  the  quality  of  solutions  contained  in  these  neighborhoods  for 
small-  to  moderate-sized  problem  instances,  and  compare  the  results  with  the  optimal  solutions  where  it 
is  possible. 

Before  proceeding  with  the  numerical  results  of  the  study,  in  the  next  section,  FINDCLIQUE,  the  al¬ 
gorithm  that  is  used  to  find  the  optimum  clique  in  the  index-graph  Q*  or  the  first  clique  in  the  a -set  or 
2a-sct  will  be  described.  The  results  from  randomly  generated  MAP  instances  for  each  of  these  two 
methods  arc  presented  next. 


4.3.1  Finding  //-Cliques  in  « -Partite  Graphs 

In  order  to  find  cliques  in  Q* ,  the  a-set,  or  the  2a-set,  the  branch-and-bound  algorithm  proposed  in  ?  is 
used.  This  algorithm,  called  FINDCLIQUE,  is  designed  to  find  all  //-cliques  contained  in  an  unweighted 
/7 -partite  graph. 

The  input  to  original  FINDCLIQUE  is  an  //-partite  graph  G  ( fj . If:  E)  with  the  adjacency  matrix 

M  =  (ntij),  and  the  output  will  be  a  list  of  all  //-cliques  contained  in  G.  Nodes  from  G  arc  copied 
into  a  set  called  compatible  nodes,  denoted  by  C .  The  set  C  is  further  divided  into  n  partitions, 

each  denoted  by  C;-  that  arc  initialized  such  that  they  contain  nodes  from  part ite  V,,  i  =  {1 . //}. 

FINDCLIQUE  also  maintains  two  other  sets,  namely,  current  clique,  denoted  by  Q  and  erased 
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nodes,  denoted  by  E.  The  set  Q  holds  a  set  of  nodes  that  are  pairwise  adjacent  and  construct  a  clique. 
The  erased  node  set,  E ,  is  furthered  partitioned  into  n  sets,  denoted  by  Ei,  that  arc  initialized  as 
empty.  At  each  step  of  the  algorithm,  Ei  will  contain  the  nodes  that  arc  not  adjacent  to  the  zth  node 
added  to  Q . 

The  branch-and-bound  tree  has  //  levels,  and  FINDCLIQUE  searches  for  // -cliques  in  the  tree  in  a  depth- 

first  fashion.  At  level  t  of  the  branch  of  bound  algorithm,  the  index  of  the  smallest  partition  in  C ,  6  = 

arg  min{  |  C,  \  ■  i  V}  will  be  detected,  and  Cg  will  be  marked  as  visited  by  including  6  into  V  { V  U  0 }, 
i 

where  V  is  the  list  of  partitions  that  have  a  node  in  O .  Then,  a  node  q  from  Cg  is  selected  at  random  and 
added  to  Q.  If  \  Q\  =  n,  an  //-clique  is  found.  Otherwise,  C  will  be  updated;  every  partition  C,  where 

i  £  V  will  be  searched  for  nodes  c/y ,  ( j  =  1 . \Ct  |)  that  are  not  adjacent  to  q,  i.e.  mq^Cij  =  0.  Any 

such  node  will  be  removed  from  Cj  and  will  be  transferred  to  Et.  Note  that  in  contrast  to  C,  nodes  in 
different  levels  of  E  will  not  necessarily  be  from  the  same  part ite  of  G .  Decision  regarding  backtracking 
is  made  after  C  is  updated.  It  is  obvious  that  in  an  //-partite  graph  the  following  will  hold: 

co(G)  <  n,  (48) 

where  o>(G)  is  the  size  of  a  maximum  clique  in  G.  In  other  words,  the  size  of  any  maximum  clique 
cannot  be  larger  than  the  number  of  partites,  in  that  the  maximum  clique  can  only  contain  at  most  1 
node  from  each  partite  of  G.  If  after  updating,  there  is  any  C,  ^  V  with  |C,-|  =  0,  adding  q,  to  Q  will 
not  result  in  a  clique  of  size  n,  since  the  condition  in  (48)  changes  into  strict  inequality.  In  such  cases, 
q  is  removed  from  Q,  nodes  from  Et  will  be  transferred  back  to  their  respective  partitions  in  C,  and 
FINDCLIQUE  will  tty  to  add  another  node  from  Cg  that  is  not  already  branched  on,  to  Q.  If  such  a  node 
does  not  exist,  the  list  of  visited  partitions  will  be  updated  (V  M\9),  and  FINDCLIQUE  backtracks 
to  the  previous  level  of  the  branch-and-bound  tree.  If  the  backtracking  condition  is  not  met  and  q  is 
promising,  FINDCLIQUE  will  go  one  level  deeper  in  the  tree,  finds  the  next  smallest  partition  in  the 
updated  C  and  tries  to  add  a  new  node  to  Q . 

When  solving  the  clique  problem  in  the  a-set  or  2a-sct,  since  the  objective  is  to  find  the  first  //-clique 
regardless  of  its  cost,  FINDCLIQUE  can  be  used  without  any  modifications,  and  the  weights  of  the  nodes 
in  CL/*)  or  Q*nn(2a)  will  be  ignored.  However,  when  the  optimal  clique  with  the  smallest  cost  in  Q* 
is  sought,  some  modifications  in  FINDCLIQUE  arc  necessary  to  enable  it  to  deal  with  weighted  graphs. 
The  simplest  way  to  adjust  FINDCLIQUE  is  to  compute  the  weight  of  the  //-cliques  as  they  are  found, 
and  report  the  clique  with  the  smallest  cost  as  the  output  of  the  algorithm.  This  is  the  method  that  is 
used  in  the  experimental  studies  whenever  the  optimal  solution  is  desired.  However,  to  obtain  a  more 
efficient  algorithm,  it  is  possible  to  calculate  the  weight  of  the  partial  clique  contained  in  Q  in  every  step 
of  the  algorithm  and  fathom  subproblems  for  which  Wq  >  Wq*,  where  Wq  and  Wq*  arc  the  cost  of 
the  partial  clique  in  O  and  the  cost  of  the  best  clique  found  so  far-  by  the  algorithm  respectively.  Further 
improvement  can  be  achieved  by  sorting  the  nodes  in  Cj,  i  —  1 ,...,//,  based  on  their  cost  coefficients, 
and  each  time  select  the  untraversed  node  with  the  smallest  node  as  the  next  node  to  be  added  to  Q 
(as  opposed  to  randomly  selecting  a  node,  which  does  not  change  the  overall  computational  time  in  the 
unweighted  graph  if  a  list  of  all  //-cliques  is  desired).  This  enables  us  to  compute  a  lower  bound  on  the 
cost  of  the  maximum  clique  that  the  nodes  in  Q  may  lead  to  as  follows: 

LBq  =  WQ  +  J2wfn-  (49) 

iyv 

where  u/-11"1  is  the  weight  of  the  node  with  the  smallest  cost  coefficient  in  C, .  Any  subproblem  with 
LBq  >  Wq*  will  be  fathomed. 
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4.3.2  Random  Linear  MAPs  of  Large  Cardinality 


To  demonstrate  the  performance  of  the  method  described  in  section  4.1,  random  MAPs  with  fixed  di¬ 
mensionality  d  =  3  and  different  values  of  cardinality  n  arc  generated.  The  cost  coefficients  <plr.. [d 
are  randomly  drawn  from  the  uniform  U  [0,  1]  distribution.  Three  sets  of  problems  arc  solved  for  this 

case:  (i)  n  =■  3 . 8  with  d  =  3,  solved  for  optimality,  and  the  first  clique  in  the  a-  and  2a-sets, 

(ii)  n  —  10,  15, ... ,  45,  with  d  =  3,  solved  for  the  first  clique  in  the  a-  and  2a-sets,  and  finally  (iii) 
n  =  50,  55, . . . ,  80,  with  d  =  3,  solved  for  the  first  clique  in  the  2a-set.  For  each  value  of  n,  25  in¬ 
stances  are  generated  and  solved  by  modified  FINDCLIQUE  for  the  optimum  clique  or  FINDCLIQUE 
whenever  the  first  clique  in  the  problem  is  desired.  Algorithm  is  terminated  if  the  computational  time 
needed  to  solve  an  instance  exceeds  1  hour. 

In  the  first  group,  (i),  instances  of  MAP  that  admit  solution  to  optimality  in  a  reasonable  time  were  solved. 
The  results  from  this  subset  arc  used  to  determine  the  applicability  of  Corollary  3.1  and  bounds  (44)  and 
(45)  for  relatively  small  values  of  n.  Table  2  summarizes  the  average  values  for  the  cost  of  the  clique  and 
computational  time  needed  for  MAPs  with  the  linear  sum  objective  function  for  the  instances  in  group  (i). 
The  first  column,  n ,  is  the  cardinality  of  the  problem.  The  columns  under  the  heading  “Exact”  contain 
the  values  related  to  the  optimal  clique  in  Q* .  The  columns  under  the  heading  “(7*in(aw)”  represent 
the  values  obtained  from  solving  the  a-set  for  the  first  clique,  and  those  under  the  heading  “t/*in(2a)” 
represent  the  values  obtained  from  solving  the  2a-sct  for  the  first  clique.  For  each  of  these  multicolumns, 
T  denotes  the  average  computational  time  in  seconds,  Z  is  the  average  cost  of  the  cliques,  |  V\  is  the  order 
of  the  graph  or  induced  subgraph  in  Q* ,  (7*in(a),  or  (/*in  (2a),  and  3  CLQ  shows  the  percentage  of  the 
problems  for  which  the  a-set  or  2a -set,  respectively,  contains  a  clique.  This  value  is  100%  for  the  exact 
method.  There  was  no  instances  in  group  (i)  for  which  the  computational  time  exceeded  1  hour. 

It  is  clear  that  using  a-sct  or  2a-sct  enables  us  to  obtain  a  high-quality  solution  in  a  much  shorter  time  by 
merely  searching  a  significantly  smaller  paid  of  the  index  graph  Q*.  Based  on  the  values  for  Z,  the  cost 
of  the  clique  found  in  a-set  or  2a-sct  arc  consistently  converging  to  that  of  the  optimal  clique  and  they 
provide  tight  upper  bounds  for  the  optimum  cost.  Additionally,  as  is  shown  in  the  \V\  column,  significant 
reduction  in  the  size  of  the  graph  can  be  obtained  if  a-set  or  2a-set  are  used. 

Table  3  contains  the  corresponding  results  for  the  case  of  a  random  MAP  with  bottleneck  objective.  In 
this  table,  W  represents  the  value  for  the  cost  of  the  optimal  clique  or  the  first  clique  found  in  a-  or 
2a-set.  Figure  4(a)  shows  how  the  cost  of  an  optimum  clique  compares  to  the  cost  of  the  clique  found 
in  a-set  and  2a-set.  Clearly,  the  cost  of  optimal  clique  approaches  0  for  both  linear  sum  and  linear 
bottleneck  MAPs.  Figure  4(b)  demonstrates  the  computational  time  for  instances  in  group  (i). 

The  advantage  of  using  a-set  over  2a-set  is  that  the  quality  of  the  detected  clique  is  expected  to  be  higher. 
On  average,  however,  a  clique  in  2a-set  is  found  in  a  shorter  time  than  in  a-set. 

The  second  group  of  problems,  (ii),  comprises  instances  that  cannot  be  solved  to  optimality  within  1 
hour.  The  range  of  n  for  this  group  is  such  that  the  first  clique  in  the  a-set  is  expected  to  be  found  within 
1  hour.  Tables  4  and  5  summarize  the  results  obtained  for  this  group.  Instances  with  n  —  45  were  the 
largest  problems  in  this  group  for  which  a-set  could  be  solved  within  1  hour.  As  it  is  expected,  the 
2a-set  can  be  solved  quickly  in  a  matter  of  seconds  where  the  equivalent  problem  for  a-set  requires  a 
significantly  longer  computational  time.  However,  the  quality  of  the  solutions  found  for  a-set  is  higher 
than  the  quality  for  solutions  in  2a-set.  Nonetheless,  using  2a-set  increases  the  odds  of  finding  a  clique, 
as  based  on  lemma  5,  2a-set  is  expected  to  contain  an  exponential  number  of  cliques.  It  is  obvious 
from  the  3  CLQ  column  that  not  all  of  the  instances  in  a-set  contain  at  least  a  clique,  whereas  100% 
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Table  2:  Comparison  of  the  computational  time  and  cost  for  the  optimum  clique  and  the  first  clique  found 
in  Q*(a)  and  Q*( 2a)  in  random  MAPs  with  linear  sum  objective  functions  for  instances  in  group  (i). 
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Table  3:  Comparison  of  the  computational  time  and  cost  for  the  optimum  clique  and  the  first  clique  found 


in  G*  (a)  and  Q*(2a)  in  random  MAPs  with  linear  bottleneck  objective  functions  for  instances  in  group 
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of  the  instances  in  2a-set  contain  one  that  can  be  found  within  1  hour.  Column  Timeout  represents  the 
percentage  of  the  problems  that  could  not  be  solved  within  the  allocated  1  hour  time  limit.  Out  of  25 
instances  solved  for  n  =  45,  only  4  (16%)  could  not  be  solved  in  1  hour.  Out  of  the  21  remaining 
instances,  20  instances  contained  a  clique,  and  only  1  did  not  have  a  clique.  The  behavior  of  the  average 
cost  values  for  the  problems  solved  in  this  group  are  depicted  in  Figure  5. 

Finally,  the  third  group,  (iii),  includes  instances  for  which  the  cardinality  of  the  problem  prevents  the 
a-set  from  being  solved  within  1  hour.  Thus,  for  this  set,  only  the  2a-sct  is  used.  The  instances  of  this 
group  were  solved  with  the  parameter  values  n  =  50,  55, . . . ,  80  and  d  =  3.  Tables  6  and  7  summarize 
the  corresponding  results.  When  the  size  of  the  problem  n  >  55,  some  instances  of  problems  become 
impossible  to  solve  within  1  hour  time  limit.  The  average  cost  for  the  instances  that  arc  solved  keeps  the 
usual  trend  and  converges  to  0  as  n  grows.  The  largest  problems  attempted  to  be  solved  in  this  group  arc 


Table  4:  Comparison  of  the  computational  time  and  cost  for  the  first  clique  found  in  G*(a)  and  G*(2a  ) 


in  random  MAPs  with  linear  sum  objective  functions  for  instances  in  group  (ii). 
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(a)  (b) 

Figure  4:  Solution  costs  (a)  and  computational  time  (b)  in  random  MAPs  with  linear  sum  and  linear 
bottleneck  objective  functions  for  instances  in  group  (i). 


Table  5:  Comparison  of  the  computational  time  and  cost  for  the  first  clique  found  in  G*(oe)  and  Q*{ 2a) 


in  random  MAPs  with  linear  bottleneck  objective  functions  for  instances  in  group  (ii). 
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MAPs  with  n  =  80.  Out  of  25  instances  of  this  size,  only  4  could  be  solved  within  1  hour.  Figure  5(a) 
the  average  values  of  solution  cost  and  computational  time  for  the  instances  of  both  linear  sum  and  linear 
bottleneck  MAPs.  Note  that  as  the  size  of  the  problem  increases,  the  reduction  in  the  size  of  problem 
achieved  from  using  a-set  or  2a-set  becomes  significantly  larger.  For  instance,  in  MAP  with  n  —  80  and 
d  —  3,  the  2a-set  has  80  x  14  nodes,  while  the  complete  index  graph  will  have  80  x  802  nodes. 


4.3.3  Random  MAPs  of  Large  Dimensionality 

The  second  set  of  problem  instances  includes  MAPs  that  are  solved  by  the  heuristic  method  explained  in 
section  4.2.  Problems  in  this  set  have  the  cardinality  n  —  2, ...  ,5  and  dimensionality  in  the  range  d  — 
2 , ...  ,dn,  where  dn  is  the  largest  value  for  d  for  which  an  MAP  with  cardinality  n  can  be  solved  within 
1  hour  using  the  heuristic  method.  For  each  pair  of  (n.  d),  25  instances  of  MAP  with  cost  coefficients 
randomly  drawn  from  the  uniform  T7  [0,  1]  distribution  are  generated.  Generated  instances  are  then  solved 
by  the  modified  FINDCLIQUE  for  the  optimal  clique  (when  possible)  and  the  optimal  costs  are  compared 
with  the  costs  obtained  from  the  heuristic  method.  The  result  of  the  heuristic  method  for  instances  with 
77  =  2  is  optimal,  and  the  heuristic  checks  all  the  2d~l  solutions  of  the  MAP.  Thus,  using  the  modified 
FINDCLIQUE  to  find  the  optimum  clique  is  not  necessary. 
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Table  6:  Computational  time  and  cost  for  the  first  clique  found  in  Q* (2a)  in  random  MAPs  with  linear 
sum  objective  functions  for  instances  in  group  (iii). 
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Table  7:  Computational  time  and  cost  for  the  first  clique  found  in  Q* (2a)  in  random  MAPs  with  linear 
bottleneck  objective  functions  for  instances  in  group  (iii). 
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Figure  6  demonstrates  the  cost  convergence  in  instances  with  n  =  2,  3,4,  5  for  both  linear  sum  and 
lineal-  bottleneck  MAPs.  Figure  6(a)  demonstrates  the  cost  convergence  in  MAPs  with  n  =  2  and 

d  =  2 . 27.  Recall  that  due  to  Remark  2,  for  cases  with  n  —  2  the  heuristic  provides  the  optimal 

solution.  The  heuristic  method  provides  high  quality  solutions  that  are  consistently  converging  to  the 
optimal  solution  for  all  cases  and  the  average  value  of  the  obtained  costs  from  the  heuristics  approaches 
0.  Memory  limitations,  as  opposed  to  computational  time,  were  the  restrictive  factor  for  solving  larger 
instances  as  the  computational  time  for  the  problems  of  this  set  never  exceeded  700  seconds. 

Figure  7  demonstrates  the  computational  time  for  the  optimal  method  as  well  as  the  heuristic  method  in 
instances  with  n  =  2,  3,  4,  5  for  both  linear  sum  and  linear  bottleneck  MAPs.  The  computational  time 
has  an  exponential  trend  as  the  number  of  solutions  for  the  MAP,  or  the  number  of  solutions  checked  by 
the  heuristic  grow  in  an  exponential  manner.  However,  the  heuristic  method  is  able  to  find  high  quality 
solutions  in  significantly  shorter  time. 


5  On  y -cliques  in  random  graphs 

A  complete  subset,  or  a  clique  in  a  simple  undirected  graph  G  =  (V.  E)  is  a  subset  Q  C  V  of  vertices 
that  are  pairwise  adjacent,  i.e.,  for  any  i,  j  e  Q  there  exists  an  edge  (/,  j)  e  E.  A  clique  Q  that  cannot 
be  increased  by  adding  new  vertices  from  V\Q  is  called  maximal,  a  clique  of  the  largest  size  is  called 
the  maximum  clique.  In  the  present  endeavor  we  are  concerned  with  the  properties  and  behavior  of  a 
certain  type  of  clique  relaxation,  namely  the  quasi-clique,  or  y-clique  Abello  et  al.  (2002). 

Definition  1  (quasi-clique).  Let  G  =  (V,  E )  be  a  simple  undirected  graph,  and  Q  C  V  be  a  subset  of 
its  vertices.  The  (induced)  subgraph  G[Q\  is  called  a  y-clique  for  a  given  fixed  y  e  (0,  1],  if  the  ratio  of 
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Figure  5:  Comparison  of  the  cost  (a)  and  computational  time  (b)  for  MAPs  with  linear  sum  and  linear 
bottleneck  objective  functions  for  group  (ii)  and  (iii). 


the  number  of  edges  in  G[Q ]  to  the  maximum  possible  number  of  edges  among  vertices  in  Q  is  at  least 


y: 


\E(Q)\>y 


Note  that  the  case  of  y  —  0  would  be  trivial,  as  any  graph  is  a  0-clique.  A  complete  graph  is  a  1 -clique, 
hence  y-clique  represents  a  density-based  relaxation  of  the  clique,  as  compared  to  degree-  and  diameter- 
based  clique  relaxations  such  as  k-plex  and  k-club  (Balasundaram  et  ah,  2011;  Alba,  1973;  Mokken, 
1979;  Wasserman  and  Faust,  1994;  Scott,  2007).  Similarly  to  maximum  cliques,  a  y-clique  with  the 
largest  number  of  vertices  is  called  the  maximum  y-clique. 

In  this  work,  we  investigate  the  asymptotical  behavior  of  y-cliques  in  large-scale  random  graphs,  and 
develop  a  compact  linear  mixed-integer  programming  formulation  for  identifying  the  largest  y-clique  in 
a  given  network. 

In  particular,  we  employ  a  popular  G  ( n ,  p)  model  of  random  graphs,  originated  by  Erdos  Erdos  (1947), 
which  denotes  a  graph  on  n  vertices,  such  that  an  edge  between  any  two  vertices  exists  with  a  probability 
p  e  (0,  1],  independently  from  other  edges.  The  G{n,  p)  model  is  related  to  the  Erdos-Renyi  G (n .  M) 
model  of  random  graphs,  in  which  graphs  with  n  vertices  and  M  edges  are  uniformly  equiprobable  with 

probability  ^  j  .  For  the  G(n ,  p)  model  yields  graph  instances  with  a  rather  “uniform”  structure,  as 
opposed  to,  say,  power-law  graphs,  it  is  sometimes  called  a  uniform  random  graph  model  Chung  et  al. 
(2001),  a  terminology  that  will  also  be  used  in  this  paper. 

Random  graphs  and  related  structures,  such  maximum  cliques  in  random  graphs,  have  been  studied 
intensively  in  last  decades  Erdos  and  Renyi  (1960);  Bollobas  and  Erdos  (1976);  Bollobas  (2001).  One 
of  the  earliest  works  on  the  asymptotical  behavior  of  maximum  clique  in  uniform  random  graphs  is 
due  to  Matula  Matula  (1970),  who  showed  that  the  maximum  clique  size  has  a  strong  peak  around 
2  In  n/ ln(  1  / p).  Grimmett  and  McDiarmid  Grimmett  and  McDiarmid  (1976)  proved  that  as  n  — >■  oc  the 
maximum  clique  size  in  a  uniform  random  graph  G(n,  p)  is  equal  to  2  In  n/  1  n (1  / p)  +  (9  (In  In  n)  with 
probability  one. 

In  Section  2,  we  present  a  generalization  of  this  result  for  the  case  of  maximum  y-cliques.  Furthermore, 
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Figure  6:  Comparison  of  the  cost  obtained  from  the  heuristic  method  with  the  optimum  cost  in  MAPs 
with  linear  sum  and  linear  bottleneck  objective  functions  with  (a)  n  —  2,  (b)  n  =  3,  (c)  n  =  4,  and  (d) 
n  —  5 


we  demonstrate  that  the  size  of  maximum  y-cliquc  in  G(n ,  p)  undergoes  a  phase  transition  when  the 
value  of  p  is  varied  in  the  vicinity  of  the  (fixed)  value  y  e  (0,  1),  manifested  in  a  sudden  and  drastic 
change  of  size  of  the  maximum  y-cliquc  in  G(n ,  p  )  relative  to  the  size  of  the  graph  itself.  The  phenomena 
of  phase  transition  in  random  structures  are  well  known  in  many  fields  of  science  and  engineering;  some 
of  the  relevant  works  include  Bollobas  et  al.  (2005,  2007);  Luczak  (1996);  Luczak  and  Luczak  (2006); 
Luczak  (1990);  Luczak  et  al.  (1994);  Frank  and  Martel  (1995);  Krishnamachari  et  al.  (2001). 

In  Section  3,  we  develop  a  linear  mixed  integer  programming  (MIP)  formulation  for  the  problem  of 
identifying  the  maximum  clique  problem  in  a  given  graph.  In  comparison  to  traditional  formulations 
existing  in  the  literature  for,  e.g.,  maximum  clique  problem,  and  involving  either  a  non-convex  quadratic 
constraints  or  a  number  of  linear  constraints  that  is  quadratic  in  the  size  of  the  graph,  our  formulation  is 
linear  and  employs  number  of  variables  and  constraints  that  is  linear  in  the  size  of  the  graph. 


5.1  Asymptotic  behavior  of  y  -cliques  in  uniform  random  graphs 

Define  N £  as  the  (random)  number  of  y-cliques  of  size  k  in  G(n .  p).  Noting  that  there  are  (k)  different 
subgraphs  of  size  k  in  this  graph,  let  ij  be  the  indicator  variable  such  that  ij  =  1  if  the  j  -th  subgraph 
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(c)  (d) 

Figure  7 :  Comparison  of  the  computational  time  in  logarithmic  scale  needed  for  the  optimal  method  and 
the  heuristic  method  in  MAPs  with  linear  sum  and  linear  bottleneck  objective  functions  with  (a)  n  =  2, 
(b)  n  =  3,  (c)  n  —  4,  and  (d)  n  —  5 


of  size  k  is  a  y-clique,  j  —  1, . . . ,  (^),  and  ij  —  0  otherwise,  which  allows  us  to  express  as 

© 

Nl  =  HIJ  ■ 

7  =  1 


(50) 


Obviously,  the  unconditional  probabilities  P {lj  =  1}  that  any  subgraph  of  size  A  is  a  y-clique  are  equal, 
whence  the  expected  number  of  y -cliques  of  size  k  in  G (n,  p)  is  given  by 

(A 

E[A^]  =  I  P{a  subgraph  of  size  k  in  G  is  a  y-clique} 


'A  f.  (& 

h-  I  \  m 


Pm(  1  -  /?)(*)“ 


(51) 


^=[y©] 


'  n  \  r 


j  -  Bin  (LKSJ: 
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where  Bin(/r;  n,  p)  is  the  c.d.f.  of  the  binomial  distribution.  As  it  will  be  seen,  the  integer  k  =■  kf  such 
that 

ElNkJ  = 1  <52> 

for  large  values  of  n  plays  a  central  role  in  the  sequel.  The  next  proposition  takes  a  first  step  in  evaluating 
ky 

Proposition  4 .  If  p  <  y,  the  integer  k  =  kf  that  satisfies  E \ Nfy  1  =  1  increases  with  n  in  such  a  way 
that  kf  =  o  {n),  n  1. 


Proof.  From  expression  (51)  it  is  evident  that  k  —  kvn  cannot  be  bounded  for  large  values  of  n,  since  in 
that  case  the  right-hand  side  of  (51)  would  be  equal  asymptotically  to  0(nk).  To  verify  that  kf  =  o(n), 
we  construct  an  upper  bound  on  the  right  hand  side  of  equation  (51).  Using  Stirling’s  approximation,  the 
binomial  coefficient  in  (51)  can  be  bounded  as 


(53) 


To  bound  the  summation  term  in  (51),  we  use  Chernoff’s  bound  for  the  tail  of  the  binomial  distribution 
Erdos  and  Spencer  (1974): 


tfVa-r) 

i  =m  \  / 


where  m  >  np.  In  our  case  n  —  (if),  m  =  y  (k)  (for  simplicity,  we  use  m  =  y( and  since  p  <  y, 
then  m  >  np:  thus,  the  requirement  on  m  is  valid.  Thus, 

© 

\pl(l  -  p)®-'  < 


(54) 


Combining  the  upper  bounds  in  (53)  and  (54),  we  have  that  if  k  =  kf  satisfies  E[/V7]  =  1,  whereby  the 
following  must  hold  for  large  enough  values  of  n : 

,1© 


ls(i 


fen 


1  -  p 


1 


Y 


t— y 


(55) 


Taking  logarithm  of  the  right  hand  side  of  the  above  inequality  and  dividing  by  k2.  we  obtain 

In  c, 


1  In  n 
k  +  k 


In  k  1  /  1 

k  +  2  V  k 


where  the  constant  c  has  the  form  c  = 


It  is  easy  to  see  from  the  inequality  for 


arithmetic  and  geometric  means  that  c  G  (0,  1)  for  0  <  p  <  y  <  1.  Then,  if  k  grows  with  n  such  that 
=  o(l),  the  above  expression  becomes  negative  for  sufficiently  large  n,  thereby  contradicting  the 
constructed  upper  bound  (55).  This  implies  that  kf  =  o(n),  which  proves  the  proposition.  ■ 
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Proposition  5 .  If  p  <  y,  the  integer  k  —  k„  that  satisfies  the  equality  E^/Vjy]  =  1,  is  given  by 

fori 


kn  ~ 


In 


GO' (SO 


—  -.  Inn  +  O(lnlnn),  n  ^  1. 


(56) 


In  establishing  Proposition  5  we  rely  on  the  following  result  due  to  McKay  McKay  (1989). 

Theorem  4  (McKay  McKay  (1989)).  Let  p  e  (0,  1)  be  fixed,  and  pv  <  k  <  v  for  some  v  >  1.  Define 
x  =  l<~(TPV ,  where  a  =  y/vp(  1  —  p).  Then 

E  -P)v_i  =  jV_1(l -PY~K  ^^expl €(v,K,p)/a},  (57) 


where 


0  <  e(y,K,p )  <  min{ y/n/ 8,  x  1}, 


a/wf  <J>(x)  nnt/  </>(;r)  are  ?/ie  cumulative  and  probability  density  functions  of  the  standard  normal  distri¬ 
bution,  respectively. 


Proof  of  Proposition  5.  Using  the  notations  of  Theorem  1,  let 

v  =  (2),  k  —  yv,  o  =  y/  p(  1  —  p)v,  x  =  —  v 


=  „l/2 

-p) 


where  we  note  that  k  ss  \yv]  >  pv  for  large  enough  v,  then  the  last  term  in  (57)  satisfies 
exp  {e(v,  k,  p)/o}  =  exp{0(v-1)}  =  1  +  0(v~l),  n  — >  00. 

From  the  fact  that  x  increases  with  n  (cf.  Proposition  4),  it  follows  that 


=  I  +  0(0  =  MEEV-W  (1  +  0(„-.)) ,  ;;  — *  oo, 

(f>{x)  x  y  —  p 


where  the  well-known  expansion 


1  -  0(f)  =  —}=e  ^/2  (f  1  +  0(t  3))  ,  t 
s/2  n 


00, 


was  used.  Invoking  Stirling’s  expansion  for  T(z), 


2  n 


T(z)  —  zze  2(l  +  0(z  1)),  z  ->■  00, 


we  obtain 


v  -  1 


jv-l  V2^V  1  -  y 


[yr  (i  -  y)l~v]  yv“1/2(i  +  0(v_1)). 
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Thus,  finally,  the  tail  of  the  binomial  distribution  in  (51)  can  be  estimated  as 


E 


i  =yv 


1  i  —  p 

Y  ~  P  V  1  ~Y 


y~ v-1/ 2 


i  -  p 
1  -  y 


l-y' 


(1  +  OCv-1)) , 


where  v  —  (2).  Consequently,  equation  (51)  can  asymptotically  be  written 


as 


1  1  —  p 


[k  V2nY~  P\  l-y\2 


-1/2  r 


Y 


1  ~  P 
1  -  y 


i-y"|(£) 


(1+0  (or1)) 


=  1. 


Taking  the  logarithm  of  both  sides  of  the  last  equality,  we  obtain 


/  1 \  1  (k\  (k\  (  ( k2  1 

(k  +  2 )  In  k  +  k  +  k  In  n  +  In  ci In  I  +1  In  c2  =  O  I  max  <  — ,  — 

2  \2/  \ 2 /  \  [  n  k 


where 


ci  = 


1  1  —  p 

2  n  y  —  P 


i-y 


c2  =  1  - 

.r 


i-p 

1  -  y 


i-y 


Taking  into  account  that 


-  In  (  )  =  Ink  —  In  \fl  +  0(k  J), 


(58) 


(59) 


equation  (59)  can  be  written  as 
,2 


—  Inc2  —  k  Ink  +  k(ln«  +  1  —  |  lnc2)  —  |  Ink  +  In  \flc\  =  O  ^max  |  — ,  —  . 

To  obtain  the  main  term  of  the  asymptotical  approximation  of  the  solution  of  the  last  equation,  let  us 
restate  it  in  the  form 


1 


k  _lnc2  +  — -  — + 


In  n  Ink  1  — ilnc2  3  Ink  lnV2ci 


2k2 


+ 


k2 


k2  1 
=  O  max  { — ,  — 
n  k 


In  view  of  the  fact  that  k  =  o(n)  due  to  Proposition  4,  the  above  expression  can  be  further  rewritten  as 

1  In  77  /  In  77  \ 

-lnC2  +  —  +  0i  —  \  =  0(l), 


whence  we  have  that 


k  = 


2  In  77 


In  c 


—  +  x(n)>  where  x(n)  —  o  (In  77 ) . 


To  determine  the  order  of  the  term  /(n),  we  restate  the  last  equation  as 


—  lnc2  +  In  77  —  Ink  =  0(1). 


Writing  down  the  expression  for  k  =  k„  in  the  form 

7  V  2  In  77 

kvn  =  ~ — rr  +  Xin), 

lnc2 

and  substituting  it  in  the  last  equation,  we  obtain  that  y(7?)  =  0(ln  In  77),  which  furnishes  the  statement 
of  the  proposition.  ■ 
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Next,  we  demonstrate  that  the  number  k„  (56),  which  solves  the  equation  E [ ]  =  1,  with  probability 
1  represents  an  upper  bound  on  the  size  of  the  maximum  quasi-clique  in  a  uniform  random  graph  when 
n  — >  oo.  For  this,  we  need  the  following  property  of  y-cliques. 

Proposition  6.  If  graph  G  —  (V,  E),  where  \  V\  =  n,  is  a  y -clique  for  some  fixed  y  e  (0,  1],  then  for 
any  s  <  n  there  exists  a  y -clique  of  size  s  in  G. 


Proof  For  y  —  1,  this  property  is  trivial.  In  the  case  of  y  e  (0,  1),  it  suffices  to  show  that  the  statement 
of  the  proposition  holds  for  s  =  n  —  1.  Since  G  —  (F,  E )  is  a  y-clique,  then 

Assume  that  there  exists  a  vertex  i  e  F  with  degG(/)  <  y(n  —  1).  Let  V\  =  V  \  i;  then  the  induced 
subgraph  G[F]  =  (F,  Ef)  is  also  a  y-clique,  since 


If  there  is  no  such  a  vertex,  i.e.,  degG(/)  >  y(n  —  1)  for  all  i  €  F,  then  let  j  =  arg  min;el/  degG(/)  be 
the  vertex  of  G  with  the  smallest  degree,  degG(y )  =  m  >  y(n  —  1).  As  before,  denote  Vj  =  V\j,  and 
observe  that  the  cardinality  of  the  set  of  edges  Ej  of  the  induced  subgraph  G[Vj ]  =  (Vj ,  Ej)  satisfies 


I  Ej 


2  J2  degG[F.](0  +  ^  d^G(Vj](i) 

“  ieVj-.(iJ)€E  i€Vj-.(i,ME 

>  -m(m  —  1)  +  —  (n  —  m  —  1  )m 


=  —  (n  —  2 )m  >  y 


thus  verifying  the  statement  of  the  proposition.  ■ 

Proposition  7.  Let  Mf  denote  the  ( random )  size  of  the  maximum  y-clique  in  a  uniform  random  graph 
G(n,p ).  If  0  <  p  <  y  <  1,  then 


r  2 

lim  sup  -  <  -  a.  s. , 

n^o o  Inn  lnl  /  p* 


(60) 


where  p* 


1  ~  P 
1  -  y 


t— y 


Proof  First,  observe  that 


P{MZ  >  k}  =  P{Nvk  >  1}  <  E[N£], 


(61) 
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where  the  equality  is  due  to  Proposition  6.  Define  a  sequence  k„  —  - Inn,  n  >  1;  then,  from 

In  1  /  p* 

expression  (51)  for  E[/V,f  ],  one  obtains  by  following  the  steps  in  Proposition  5  that  for  sufficiently  large 
values  of  n 


P {Ml  >  kn}  < 


r  ikn  +  1) 


(p*)c 


kn') 


(62) 


From  the  definition  of  k„  it  follows  that  the  term 

nkn~l(p*)(k2  ) 


is  bounded  for  large  enough  n,  whence  the  sought  probability  can  be  subsequently  bounded  as 


P {M%  >kn}< 


n 


< 


ne 


T(kn  +  l)~  knk" 


Again  recalling  the  definition  of  k„ ,  we  note  that 

2  n ekn  _  exp  {kn(l  +  §  In  1/p*)} 
kl  r  r 


0,  n  oc, 


(63) 


which  implies  that  P{A/,f  >  kn  j  =  o(n  2)  for  n  1,  whereby  the  upper  bound  (60)  on  the  size  of  the 
maximum  y -clique  holds  almost  surely  by  virtue  of  the  Borel-Cantelli  lemma.  ■ 


The  next  corollary  shows  that  in  sufficiently  large  random  graphs  G(n,  p),  the  size  of  the  maximum 
y-cliquc  is  almost  surely  above  a  certain  value  of  the  order  of  Inn.  It  uses  a  well  known  fact,  established 
by  Grimmett  and  McDiarmid  Grimmett  and  McDiarmid  (1976),  that  the  size  of  the  maximum  clique  in 
a  uniform  random  graph  G(n,  p)  converges  almost  surely  to  2  In  n /  In  /A1 .  Note  that  Mt]  represents  the 
size  of  the  maximum  clique  (y  =  1)  in  a  uniform  random  graph  G(n .  p).  Observe  also  that,  according 
to  (56), 

..  2  In  n 

lim  H  =  - - - 1-  6>(lnlnn)i 

y->t-o  In  1/p 

which  corresponds  to  the  well-known  expression  for  the  size  of  the  maximum  clique  in  uniform  random 
graphs  Bollobas  and  Erdos  (1976);  Grimmett  and  McDiarmid  (1976).  This  allows  us  to  define  A/  as  the 
limiting  value  of  above. 

Corollary  4.1.  If  0  <  p  <  y  <  1,  then  the  size  of  the  maximum  y -clique  in  a  uniform  random 
graph  Gin,  p)  satisfies 


2 

In  1/p 


<  lim  inf  — -  <  lim  sup  -  < 

n-+o o  In  n  n->-oo  In  n 


a.  s. 


(64) 


where  kj,  is  given  by  (56). 


Proof  This  follows  immediately  from  Proposition  7,  and  the  observation  that,  for  any  y  <  1,  the  size  of 
the  maximum  y -clique  in  G(/i,  p)  always  greater  than  the  size  of  the  maximum  clique  in  the  same  graph, 


i.e., 

P{Ml<Ml)=L  y  <  l. 
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Also  note  that  the  relations  0  <  p  <  y  <  1  imply  that 


y  i-py-y 

P  i  -  y) 


>  1  >  y, 


from  which  one  infers  the  inequality 


1  -P 
1  -  y 


t -v 


>  P- 


verifying  that  inequality  for  the  lower  and  upper  bounds  on  M„  /  In  n  in  (64)  always  holds,  given  the 
above  assumptions  on  the  values  of  p  and  y.  ■ 

Remark  3.  From  Fatou’s  lemma  it  follows  that  bounds  (64)  on  the  maximum  y-clique  size  A/,/,  which 
hold  with  probability  1,  also  hold  for  the  average  maximum  y-cliquc  size,  E [M%]: 


2 

in  i  ip 


<  liminf 

^OO 


E [M%] 
Inn 


<  hm  sup 
n^-o o 


(65) 


In  such  a  way,  we  have  established  that  for  any  fixed  y  >  p  the  asymptotic  size  of  the  maximum  y-clique 
is  of  the  order  of  In  n .  Intuitively,  when  y  <  P ,  the  entire  graph  G(n,  p)  becomes  a  y-clique,  thus  the 
size  of  the  maximum  y-clique  has  the  order  of  n.  Therefore,  the  natural  question  arising  here  is  what 
happens  when  y  is  fixed  and  p  approaches  y.  We  show  that  there  is  a  first  order  phase  transition  in  the 
asymptotic  behavior  of  the  order  of  magnitude  of  the  maximum  y-clique  in  the  point  y  =  p. 

Proposition  8.  If  M„  is  the  size  of  the  maximum  y-clique  in  a  uniform  random  graph  G(n ,  p)for  some 
fixed  y  €  (0,  1),  then  with  high  probability  (w.h.p.) 

lim  lim  -  =  0,  (66a) 

p/y  n-*o O  n 

but 

Mf 

lim  lim  -  =  1.  (66b) 

p\yll-+00  U 


Proof.  The  first  limiting  case  follows  from  Proposition  7,  since  we  proved  that  for  any  fixed  y  >  p  with 
probability  1 


lim 


Mi 


=  0. 


n  — >oo  n 


To  prove  the  equality  in  (66b),  let  Xij  be  a  Bernoulli  random  variable  which  is  equal  to  1  if  there  exists 
an  edge  (i.  j )  in  the  uniform  random  graph  G(n,  p).  Then, 


r  /  x  i 

E  Xu 

Tl 

J? 

a  ^ 

II 

3 

II 
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From  the  weak  law  of  large  numbers  it  follows  that  for  any  fixed  s  >  0, 


' 

E  Xij 

0  J) 

G)  P 

<  £ 

, 

1 ,  n  —*■  oo  . 


Letting  s  —  p  —  y,  we  obtain 


E  Xu 

E  x^ 

(hj)  . 

>  —  p  . 

(hi) 

n  -r  r« 

(2) 

>  —  r  < 
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E  Xij 
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p 
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>  P 


Therefore, 

P {M%  =  n }  — »■  1,  n  — >  oo, 
which  ends  the  proof  of  the  proposition. 


1,  n  — >  oo. 


Remark  4.  The  phase  transition,  a  phenomenon  of  a  drastic  change  in  some  property  of  a  random 
structure  over  a  small  change  in  the  structure’s  parameters,  is  well  known  in  the  literature.  With  respect 
to  random  graphs,  the  limiting  probability  of  a  graph’s  property  changing  from  0  to  1  or  vice  versa  is 
well  known  for  monotone  and  first  order  graph  properties  Alon  and  Spencer  (2000).  Recall  that  property 
Q  is  monotone  increasing  (respectively,  decreasing )  if  from  A  C  B  (resp.,  B  C  A)  and  A  e  Q  it 
follows  that  B  €  Q.  The  first  order  graph  properties  arc  ones  that  can  be  finitely  described  in  a  first  order 
language,  i.e.,  language  consisting  of  variables  that  represent  graph  vertices,  equality  (=)  and  adjacency 
(~)  relations.  Boolean  symbols  V,  A,  and  the  universal  and  existential  quantifications  V,  3.  Note  that 
first  order  properties  are  not  necessarily  monotone  and  vice  versa;  for  instance,  the  increasing  property 
“graph  is  connected”  cannot  be  expressed  in  first  order  language  Janson  et  al.  (2000).  Then,  limiting 
relations  similar  to  (66)  that  concern  random  graphs  with  first  order  properties  Q  arc  known  as  zero-one 
laws  Alon  and  Spencer  (2000);  Janson  et  al.  (2000): 

lim  P {G(n,  p)  has  Q)  =  0  or  1, 

n^-o o 

where  the  probability  is  monotone  if  Q  is  monotone. 

In  this  context,  it  is  worth  noting  that  the  property  that  “graph  is  a  y-clique”  in  neither  monotone,  nor 
first  order  property,  hence  the  phase  transition  in  the  relative  size  of  y-clique  in  uniform  random  graphs 
(66)  may  not  be  obtained  directly  from  the  general  zero-one  laws  relations. 


5.2  Linear  mixed-integer  formulations  of  the  maximum  y-clique  problem 

In  this  section  we  develop  linear  mixed-integer  formulations  of  the  maximum  y-clique  problem.  To  the 
best  of  our  knowledge,  the  formulation  that  we  present  at  the  end  of  this  section  is  the  most  compact 
lineal-  formulation,  of  the  size  of  0(/i),  where  n  is  the  number  of  vertices  in  the  graph. 

Consider  a  graph  G  =  (V,  E)  with  n  vertices  and  an  adjacency  matrix  A,  and  suppose  we  select  some 
subgraph  Gs  of  G.  In  order  to  verify  whether  Gs  is  a  y-clique,  we  introduce  the  binary  vector  of  variables 
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x  €  {0,  1}”,  where  x,-  =  1  if  vertex  i  belongs  to  Gs,  and  x,-  =  0  otherwise.  The  subgraph  Gs  is  a  y-clique 
if  the  cardinality  of  its  set  of  edges  is  at  least 


Y 


1  n  /  n  \  \  (  n  n 

2  Y  £  Xi  (  £  Xi  ~  1 )  =  2y\  S  XJXl  ~  Y Xi 

1  =  1  \i  =  l  /  \i,j= l  i  =  l 

1  /  n  n  n  \ 

v'  L  Av-V'  +  £  v<2  -  I> 

Viw1  i  =  1  ) 


-y  xi xi > 


ij  =  1 


where  the  last  equality  is  due  to  x;2  =  x; .  The  number  of  edges  in  the  subgraph  Gs  can  be  calculated  as 


-x‘ Ax  —  Y,  aijXjXj. 

i,j  =  i 
i¥=j 

Therefore,  the  problem  of  finding  the  maximum  y-clique  in  the  graph  G  can  be  formulated  as  follows: 


max 

xe{o,r}" 

n 

Yx> 

i  =  1 

n 

n 

s.  t. 

Y  aijxjxi  >  Y 

Y  xJxi 

i,j  =  1 

ij  =  1 

i#7 

*#7 

(67) 


This  is  a  0-1  integer  programming  (IP)  problem  with  a  linear  objective  and  a  nonconvex  quadratic  con¬ 
straint.  A  linearization  of  this  problem  can  be  performed  at  the  expense  of  introducing  additional  vari¬ 
ables  and  constraints.  To  this  end,  define  Wij  =  Xi  Xj  for  every  pair  of  nodes  (/,  j  ).  Note  that  only 
n(n  —  l)/2  —  n  of  such  variables  arc  required  since  w,j  =  u;/;.  Also,  the  constraint  =  x*xy  is 
equivalent  to 

Wij  <  Xi, 

Wij  <  Xj, 

Wij  >  Xi  +  Xj  -  1. 

Now,  we  can  reformulate  (67)  as  a  linear-  problem 


n 

max  ^  x,- 
i  =  r 

n 

s.  t.  Y,aijWij>y 

i,j  =  1 
'¥=j 

(68a) 

n 

Y  WV’ 

1.7  =  1 
*#7 

(68b) 

Wij  <  Xi, 

0  7 

=  1,.. 

,,n. 

(68c) 

K1 

VI 

3 

0  7 

=  1,.. 

. ,  n , 

(68d) 

Wij  >  +  X/  - 

1, 

0  7 

=  1,.. 

.  ,n, 

(68e) 

X$ ,  Wij  G  {0,  1}, 

0  7 

=  1,.. 

. ,  n . 

(68f) 
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This  problem  is  a  linear  0-1  problem  with  n(n  —  l)/2  variables  and  | n(n  —  1)  +  1  constraints.  We  can 
further  rewrite  it  in  a  more  compact  form  using  the  fact  that  jtj  =  0  implies  w;z/  =  0,  j  =  /  +  1 
Thus,  instead  of  (68c)  and  (68d)  we  may  write 

n 

Wij  <  nxi,  i  =  l, ...  ,n  —  l, 

j=i  + 1 

which  reduces  the  number  of  constraints  to  (n  —  1)  +  n. 

In  what  follows,  we  present  an  improved  mixed-integer  linear  formulation  of  (67)  with  only  0(n)  vari¬ 
ables  and  constraints.  Recall  that  originally  we  had  only  one  constraint, 

n 

X)  (a O'  “  Y)xJxi  >  0, 

i,y=i 

which  we  may  rewrite  as 

n  /  n  \ 

yxz  +  -  y)v7  >  0. 

t=l  V  7=1  / 

Let  us  define  the  variables  Wi  sR,i  =  1 . n,  as  follows 

(  H  \ 

w n  =  x,-  yxi  +  -  y)xj  I  ,  i  =  l,...,n. 

Next,  observe  that  each  of  the  quadratic  equalities  above  is  equivalent  to  four  linear  inequalities 

W{  <  nxi, 

Wi  >  —nxi , 
n 

Wi  >  yxz  +  y^(fl;y  _  y)^7  -  (1  -  Xj)n, 

7  =  1 
n 

Wi  <  yx,-  +  T>;/  -  y)Xj  +  (1  -  Xi)n, 

7  =  1 

X;  e  {0, 1},  Wi  e  R. 

Therefore,  the  problem  of  finding  a  maximum  y-clique  can  be  represented  as  the  following  mixed  integer 
lineal-  programming  problem  with  2 n  variables  (n  binai'y  variables  and  n  continuous  variables)  and  An  + 1 
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constraints: 


max  ^  Xi 
i= t 
n 

s.  t.  ^  Wj  >  0, 
i  =  1 


iv i  <  nxi ,  i  —  1 

Wi  >  —nxi,  i  —  1 

n 

Wi  >  yx,  +  -  y)x7  -  (1  -  x;)zz,  i  =  1 

7  =  1 
n 

Wi  <  yxi  +  y^(fl;7  -  y)x7  +  (1  -  Xi)n,  i  -  1 
7  =  1 

X;  e  {0, 1},  Wi  eR,  z  =  1 


,  Z2, 


,  ZZ, 


,  Z2, 


,  ZZ, 


,  ZJ. 


(69) 


5.3  Computational  experiments 

In  this  section  we  illustrate  the  theoretical  developments  presented  above  using  numerical  computations 
of  the  maximum  y-clique  in  randomly  generated  graphs.  In  particular',  the  numerical  studies  arc  aimed 
at  elucidating  the  following  two  aspects:  (i)  whether  the  asymptotic  bounds  (64),  (65)  on  the  size  of  the 
maximum  y-cliquc  and  its  mean  value  hold  for  relatively  small  values  of  zz,  and  (ii)  the  manifestation 
of  the  phase  transition  effect  in  the  order  of  magnitude  of  the  maximum  y-clique  in  finite-size  random 
graphs  when  p  approaches  y. 

According  to  Remark  3,  in  large  enough  random  graphs  G(n ,  p)  the  average  size  E[M„  ]  of  the  maximum 
y-clique  belongs  to  the  interval 


\mn-mn\ 


2  In  zj  2  In  zz 

Hvr  >^mn_ 


(70) 


provided  that  p  <  y.  Therefore,  it  was  of  interest  to  determine  the  applicability  of  the  above  bounds  for 
relatively  small  values  of  zz. 

To  this  end,  in  the  first  set  of  computational  experiments  we  generated  a  number  of  instances  of  uniform 
random  graphs  G(n,p)  with  zz  =  100  and  p  ranging  from  0.05  to  0.15;  namely,  we  generated  100 
instances  of  G(100,  p)  for  every  p.  Then,  we  employed  the  MIP  formulation  (69)  to  find  the  maximum 
y-cliques  in  the  generated  graphs  for  values  y  =  0.9  and  y  =  0.85.  Such  a  choice  of  parameters 
is  justified  by  relatively  better  numerical  tractability  of  the  MIP  problem  (69)  for  sparse  graphs.  We 
used  FICO™  Xpress  Optimization  Suite  7.1  Xpress  to  solve  the  resulting  instances  of  problem  (69). 
The  resulting  average  values  of  as  well  as  the  minimum  and  maximum  values  of  M^00  over  100 

instances  for  each  p,  are  reported  in  Table  8. 

In  the  second  set  of  computational  experiments,  we  analyzed  the  behavior  of  the  relative  size  of  the  max¬ 
imum  y-clique  for  a  fixed  y  and  different  values  of  p  and  n.  We  used  y  =  0.70  and  two  sequences  of  p 
values:  {p\ _ _  />ioo}  =  {0.006,  0.012 . 0.600}  and  {/?ioi, _ P200}  =  {0.601, 0.602, . . . ,  0.700}. 
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Note  the  increased  “density”  of  the  second  sequence,  which  allows  for  a  more  thorough  investigation  of 
the  maximum  y-clique  size  when  the  value  of  p  becomes  close  to  y.  For  each  value  of  n  —  500,  1000, 
5000,  1000,  20000  and  pi  from  the  sequence  defined  above,  we  generated  instances  of  uniform  random 
graphs  G(n ,  p, ),  and  determined  the  maximum  y-clique  size  M% .  Since  the  MIP  formulation  (69)  be¬ 
comes  computationally  intractable  for  large  dense  graphs,  we  used  maximum  y-clique  GRASP  heuristics 
due  to  Abello  et  al.  (2002),  which  were  reported  to  perform  quite  well  in  massive  graphs.  Figure  8  reports 
the  results  of  the  described  computational  studies,  and  Figure  9  presents  the  results  of  similar  studies  for 
y  =  0.5. 

Recall  that  it  was  shown  in  Section  2  that  with  probability  that  approaches  to  1  one  has 


lim 


n— >oo  n 


0,  p  <  y, 
1,  p  >  y, 


or,  in  other  words,  the  relative  size  of  the  maximum  y-clique  in  uniform  random  graphs  as  a  function 
of  the  density  p  of  the  graph  with  high  probability  represents  a  step  function  as  n  — >•  oo.  The  results 
presented  in  Figures  8  and  9  illustrate  the  convergence  of  the  relative  size  of  maximum  clique  to  the 
aforementioned  limits  as  n  increases. 


Figure  8:  Relative  size  of  the  maximum  y-cliques  in  the  uniform  random  graphs  for  y  =  0.7,  n  =  500, 
1000,  5000,  1000,  20000,  and  0  <  p  <  0.8. 


6  A  bit-parallel  algorithm  for  finding  k -cliques  in  a  k -partite  graph 

Given  an  (undirected)  graph  G  —  (F,  E),  where  V  is  set  of  nodes  and  E  is  the  set  of  arcs,  a  clique  in  G 
is  defined  as  a  complete  subset  of  G,  i.e.,  a  set  of  nodes  in  V  that  are  pairwise  adjacent.  A  clique  of  size 
k  is  called  k -clique ;  the  largest  clique  in  a  graph  is  called  the  maximum  clique  and  its  size  is  denoted  by 
<w(G).  Note  that  G  may  contain  several  cliques  of  size  o)(G).  Closely  related  to  the  concept  of  a  clique  is 
that  of  an  independent  set  of  G,  defined  as  an  induced  subgraph  of  V  whose  nodes  arc  pairwise  disjoint. 


43 


.  n=500 

-  n=1000 

- n=5000 

n=10000 
- n=20000 


Figure  9:  Relative  size  of  the  maximum  y-cliqucs  in  the  uniform  random  graphs  for  y  =  0.5,  n  =  500, 
1000,  5000,  1000,  20000,  and  0  <  p  <  0.6. 


Table  8:  Average  (E  minimum  (M„),  and  maximum  {Mvn)  values  of  the  maximum  y -clique  size, 

computed  over  100  instances  of  uniform  random  graphs  G(n,  p )  for  n  —  100  and  y  —  0.85,  0.90.  The 
lower  and  upper  bounds  m\,  n„  are  given  by  (70). 
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The  Maximum  Clique  Problem  (MCP)  consists  in  finding  the  largest  clique  in  a  graph,  and  is  of  fun¬ 
damental  importance  in  discrete  mathematics,  computer  science,  operations  research,  and  related  fields 
Bomze  et  al.  (1999).  In  many  applications  it  is  of  interest  to  identify  all  maximum  cliques  in  a  graph. 
This  problem  is  known  as  the  Maximum  Clique  Enumeration  Problem  (MCEP).  In  the  present  work,  we 
consider  a  special  case  of  the  MCEP,  concerned  with  finding  all  A -cliques  in  a  A -partite  graph.  A  graph 
G  —  ( V ,  E)  is  called  A -partite  if  the  set  of  nodes  V  can  be  partitioned  into  A  independent  sets,  or  partites 
Vr,  r  =  1 . k: 

k 

V  —  |J  Vr,  Vr  fl  Vs  =  0,  r  ^  s,  such  that  for  all  i,j  e  Vr  :  (i .  j )  £  E.  (71) 

r  =  1 

Clearly,  one  has  that  a>(G)  <  k  in  a  A: -partite  graph  G,  since  the  maximum  clique  cannot  contain  more 
than  one  node  from  each  independent  set  Vr.  The  problem  of  finding  A: -cliques  in  A: -partite  graphs  has 
found  applications  in  textile  industry  Grunert  et  al.  (2002),  data  mining  and  clustering  Peters  (2005),  and 
identification  of  protein  structures  Liu  and  Chen  (2009).  This  problem  is  not  necessarily  equivalent  to 
MCEP  since  it  does  not  account  for  maximum  cliques  with  co(G)  <  k. 

Grunert  et  al  Grunert  et  al.  (2002)  proposed  branch-and-bound  algorithm  FINDCLIQUE  for  the  problem 
of  finding  all  A: -cliques  in  A -part ite  graphs,  which  takes  as  an  input  a  graph  G  =  (V,  E),  where  V 
satisfies  (71),  and  produces  the  set  Q  of  A' -cliques  contained  in  G  as  an  output.  FINDCLIQUE  is  a 
recursive  method,  such  that  level  I  of  recursion  corresponds  to  the  level  t  of  branch-and-bound  tree, 
which  in  turn,  is  associated  with  the  t- th  partite  that  is  branched  on  in  V .  Starting  at  the  root  (t  =  0)  of 
the  branch-and-bound  tree  with  a  partial  solution  S  =  0,  at  each  step  of  branch-and-bound  procedure  a 
node  is  added  to  or  removed  from  S  until  S  amounts  to  a  A -clique  in  G,  i.e.,  |Sj  =  k,  or  it  is  verified 
that  G  contains  no  A:-cliques,  co(G)  <  k. 

Let  B  =  {1, . . . ,  k)  be  the  index  set  of  partites  in  G,  V  =  U beB  Vb-  and  B$  denote  the  set  of  partites 
that  have  a  node  in  S : 

Bs  =  {b  €  B  |  Vb  n  S  +  0}. 

Given  a  partial  solution  S,  a  node  is  called  compatible  if  it  is  adjacent  to  all  the  nodes  in  5;  the  set  of 
compatible  nodes  w.r.t.  S  is  denoted  by  Cs'- 

Cs  =  !>'  e  V  |  (i,j)  eE  Vj  e  Sj. 

The  set  Cs  is  further  partitioned  into  subsets  containing  nodes  from  the  same  partite: 

Cs  =  1J  Cs,b, 

beBs 


where  B s  —  B  \  Bs,  and  Cs,b  Q  Vb  is  given  by 

CS,b  =  |J<E*nA^))- 

s€S 

with  N(s)  being  the  set  of  nodes  adjacent  to  node  s. 

At  the  root  node  of  the  branch-and-bound  tree  ( t  =  0),  one  has  S  =  0,  B  =  B  s  =  {I ,  A},  Bs  =  0, 
and  Cs,b  —  Vb  for  all  b  e  B .  At  a  level  t  of  the  branch-and-bound  tree,  bt  6  B s  is  selected  as  the 
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partition  to  branch  on.  In  order  to  achieve  the  greatest  reduction  in  the  size  of  the  branch-and-bound  tree 
when  pruning,  bt  is  selected  as  the  partition  with  the  smallest  number  of  nodes: 

bt  €  argmin{|C^|  |  b  e  Bs}.  (72) 

b 

As  long  as  there  is  a  node  nt  e  Cg.i,,  that  is  not  traversed,  the  search  process  is  restarted  from  this  point 
with  S  :=  S  U  {nt}  as  the  new  partial  solution.  To  this  end,  the  set  Cs  of  compatible  nodes  is  updated 
with  respect  to  S  U  {nt}\ 

Cs,b  ■=  Cs,b  n  N(nt)  for  all  b  e  Bs.  (73) 

Maintaining  the  sets  Cs,b  of  nodes  compatible  with  the  current  partial  solution  S  is  a  key  aspect  of  the 
algorithm,  thus  for  backtracking  purposes  the  nodes  that  are  removed  from  C during  (73)  arc  added 
to  the  set  C  =  (Jf=i  which  is  similarly  partitioned  into  k  levels  C t,  each  level  corresponding  to 
level  t  of  the  branch-and-bound  tree.  In  other  words,  C  t  contains  the  nodes  in  Csj,  that  are  not  adjacent 
to  node  n  t : 

C t  —  {i  s  CSib  I  £  E,  b  €  Bs}. 

Obviously,  after  this  step,  Csj,,  =  0.  A  subproblem  with  a  partial  solution  S  is  promising  if  all  of  the 
partitions  in  Cs  that  do  not  share  a  node  in  the  partial  solution  arc  nonempty: 

\Cs,b\  >  0  for  all  b  e  Bs ,  b  ^  bt.  (74) 

Let  P  be  the  number  of  partitions  Csj,  —  Cs  that  contain  at  least  one  node;  then,  an  upper  bound  on  the 
size  of  the  largest  clique  containing  S  is  given  by  |Sj  +  P .  If  |Sj  +  P  =  k.  the  current  subproblem  is 
feasible,  meaning  S  may  be  part  of  a  A: -clique.  For  a  feasible  subproblem,  the  algorithm  traverses  deeper 
into  the  branch-and-bound  tree,  t  :=  t  +  1,  and  a  new  subproblem  is  created. 

Accordingly,  a  subproblem  with  partial  solution  S  is  pruned  if 

\S\  +  P<k ,  (75) 

i.e.,  there  exists  no  clique  of  size  k  that  contains  S.  For  a  nonpromising  subproblem,  set  Cs.b,  is  restored 
by  moving  the  nodes  in  C t  back  to  Cs.  Cs  :=  Cs  U  C t.  The  last  operation  implicitly  requires  that  the 
nodes  from  C  t  arc  put  back  into  the  partitions  of  Cs  that  they  were  removed  from: 

Cs,n{v)  • —  Cs,n(v)  hi  v  for  all  v  £  Ct ,  (76) 

where  n(i)  is  the  index  of  the  partite  that  node  i  belongs  to:  i  e  Vnjy.  moreover,  the  relative  orders 
of  nodes  in  the  partites  V},  should  be  preserved  in  Cs.b.  given  that  the  nodes  in  G  arc  assumed  to  be 
ordered/numbered. 

The  scarcli  process  is  then  restarted,  provided  that  there  exists  a  node  in  partition  Csj,,  that  is  not 
traversed.  If  there  is  no  such  node,  FINDCLIQUE  returns  to  the  previous  level  t  —  1  of  the  branch-and- 
bound  tree. 

6.1  A  bitwise  algorithm  for  finding  /v-cliques  in  a  ^-partite  graph 

In  this  section,  we  present  an  algorithm,  referred  to  as  BitCLQ,  for  the  k-clique  enumeration  problem  in  a 
k -partite  graph,  which  improves  upon  the  FINDCLIQUE  algorithm  of  Grunert  et  al  Grunert  et  al.  (2002) 
by  introducing  bitset  data  structures  and  utilizing  bit  parallelism  for  updating  the  set  of  compatible  nodes 
and  improving  backtracking. 
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6.1.1  Bitsets 


Bitsets  arc  essentially  binary  vectors,  or  sequences  of  bits,  and  as  such  can  be  utilized  efficiently  in  com¬ 
puter  codes.  Particularly,  bitsets  are  useful  for  storing  adjacency  matrices  of  graphs,  or  specific  subsets  of 
ordered  sets.  For  example,  in  a  graph  on  six  nodes  {tq , . . . ,  t>6}  =  V,  a  clique  with  nodes  vi,V2,  V3,  V5 
can  be  represented  by  a  bitset  {111010},  where  each  bit  corresponds  uniquely  to  a  node  in  the  graph,  with 
the  significant  bits  (i.e.,  bits  equal  to  1)  indicating  the  nodes  in  the  clique.  Bit  parallelism  is  a  form  of 
parallel  computing  that  achieves  computational  improvements  by  representing  the  problem  data  in  bitsets 
of  size  R .  where  R  is  the  machine  word  size  (e.g.,  32  or  64),  such  that  they  can  be  processed  together 
within  a  single  processor  instruction.  Bit  parallelism  has  been  successfully  used  in  many  computational 
algorithms,  particularly  for  string  matching  (Grabowski  and  Fredriksson,  2008;  Hyyro,  2005;  Hyyro  and 
Navarro,  2004).  Recently,  bit  parallelism  has  been  employed  for  solving  hard  combinatorial  problems, 
such  as  SAT  Segundo  et  al.  (2008)  and  the  Maximum  Clique  Problem  San  Segundo  et  al.  (201 1). 

In  the  present  work,  bit  parallelism  is  used  to  improve  the  computational  procedure  for  updating  the  set 
of  compatible  nodes  in  (73),  and,  moreover,  to  achieve  faster  backtracking  by  eliminating  the  need  for 
set  C .  In  addition,  use  of  bitsets  allows  for  improvements  in  memory  storage  efficiency  for  problem  data 
structures,  such  as  the  set  of  compatible  nodes  and  the  adjacency  matrix  of  the  graph. 

Of  particular  significance  in  the  context  of  the  present  work  is  the  operation  of  indexing  the  first 
significant  bit  in  a  bitset,  also  known  as  the  forward  bit  scanning.  One  of  the  techniques  for  this 
purpose  relies  on  use  of  the  De  Bruijn  sequence  with  a  perfect  hash  table  Leiserson  et  al.  (1998, 
http://supertech.csail.mit.edu/papers/debruijn.ps).  The  value  to  be  looked  up  in  the  hash  table  is  given  by 
Hr  below: 

Hr  \=  (x  A  -x)  D  »  (R-  log2  R),  (77) 

where  x  is  the  bitset  for  which  the  first  significant  bit  has  to  be  indexed,  D  is  an  instance  of  De  Bruijn 
sequence,  R  is  the  machine  word  size,  and  stands  for  the  binary  shift  right  operator.  Hr  is  effective 
for  bitsets  of  maximum  size  equal  to  R.  For  larger  bitsets,  special  containers  need  to  be  devised.  The 
hash  table  required  to  look  up  the  value  of  Hr  is  created  based  on  the  particular  De  Bruijn  sequence 
used  in  (77). 

Note  that  in  (77)  multiplication  is  performed  modulo  R  and  only  the  last  log2  R  bits  of  the  result  will  be 
retained.  More  details  on  forward  bit  scanning  and  the  specification  of  the  De  Bruijn  sequence  used  in 
(77)  can  be  found  in  Leiserson  et  al.  (1998,  http://supertech.csail.mit.edu/papers/debruijn.ps). 

6.1.2  BitCLQ 

Below  we  present  a  modification  of  FINDCLIQUE,  which  we  refer  to  as  BitCLQ,  that  uses  bitset  data 
structures  and  bit  parallelism  for  keeping  track  of  the  nodes  in  G  that  are  compatible  to  the  current  partial 
solution  S,  while  simultaneously  reducing  the  computational  cost  of  backtracking. 

To  this  end,  we  introduce  a  set  M  consisting  of  k  levels,  Z\, ....  Zj(.  Each  of  these  k  levels  will  be  used 
to  represent  the  compatible  nodes  to  the  partial  solution  S  at  the  7-th  level  of  the  branch-and-bound  tree, 
where  1  <  t  <  k.  Every  level  in  Z  is  further  partitioned  into  k  sets,  each  corresponding  to  a  partite  \fi 
in  G: 

—  [_J  Zt,b,  t  —  l, ...  ,k. 

beB 

The  sets  Zt  j,  arc  represented  by  bitsets  of  size  1 1/,  | .  Let  Zt  j,j  be  the  i-th  bit  in  Zt  j,  corresponding  to 
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the  z'-th  node  in  Vb,  such  that  Ztbi  =  1  if  the  z'-th  node  in  Vb  is  compatible  with  all  the  nodes  in  the 
partial  solution  S  at  the  t- th  level  of  the  branch-and-bound  tree  in  BitCLQ: 

_  (  1,  if  (/,;)  e  E  for  all  j  e  St; 
t,b’1  ~  )  0,  otherwise. 

Clearly,  each  level  Zr  of  Z  is  an  ordered  set  of  combination  of  bitsets  with  the  total  size  \V\.  Further, 
the  adjacency  matrix  M  of  graph  G  is  stored  in  the  bitset  form,  with  the  convention  that  the  i  -th  row 
(column)  corresponds  to  the  /'-th  bit  in  Zf,  t  —  1 , ,k. 

BitCLQ  is  initialized  by  setting  t  :=  0,  S  :=  0,  B  =  Bs  :=  {1 . k},  and  Q  :  =  0,  where  Q  is  the 

set  of  all  /c -cliques  in  G.  Note  that  since  at  the  beginning  all  the  nodes  in  G  can  be  added  to  S  to  extend 
its  size,  all  the  bits  in  Z\  are  significant: 

Zhbj  =  1  for  all  b  e~B(St),  i  €  Vb. 

At  level  t  of  the  branch-and-bound  tree,  the  partition  bt  to  branch  on  is  selected  as 

bt  e  argmin{|Z^|  |  b  e  Bs),  (78) 

b 

where  \Zt%b\  is  defined  as  the  number  of  significant  bits  in  the  bitset  Zt  b.  The  forward  bit  scanning 
method  discussed  in  Section  6.1.1  is  used  to  identify  node  nt  e  Vbt  that  has  not  been  traversed  and  thus 
can  be  added  to  the  partial  solution.  As  long  as  such  a  node  exists  in  Vbj,  the  search  process  is  restarted 
with  S  :=  S  U  {nt}  as  the  partial  solution,  and  the  corresponding  bit  in  Zt  bf  is  set  to  0. 

Utilizing  bitsets  also  facilitates  the  process  of  updating  the  compatible  nodes:  when  nt  is  added  to  partial 
solution,  Zt+ 1  is  created  by  performing  a  logical  AND  operation  with  Zt  and  the  row  Mint)  of  the 
adjacency  matrix  corresponding  to  the  node  nt  as  operands: 

Zt+i  =  Zt  A  M{nt).  (79) 

Similarly  to  FINDCLIQUE,  let  P  denote  the  number  of  partitions  Zt  b  with  \Zt  b\  >  0  at  level  the  t 
of  the  branch-and-bound  free.  If  |Sj  +  P  =  k,  the  current  partial  solution  is  promising,  so  that  a  new 
subproblem  is  created,  and  BitCLQ  proceeds  one  level  deeper  into  the  branch-and-bound  tree,  t  :=  t  +  1. 
If  the  partial  solution  is  not  promising,  the  method  presented  in  Section  6.1.1  is  used  to  select  nodes  in 
Vbt  that  have  not  been  traversed.  If  such  a  node  is  found,  the  search  process  is  restarted,  otherwise 
backtracking  is  performed  by  simply  updating  t  :=  t  —  1.  Note  that  due  to  the  special  structure  of 
Z,  BitCLQ  does  not  need  to  restore  the  set  of  compatible  nodes  during  backtracking,  in  contrast  to  the 
update  procedure  (76)  for  the  set  Cs  that  is  performed  in  FINDCLIQUE. 


6.1.3  Example 

As  an  illustration,  consider  the  3-partite  graph  that  is  shown  along  with  its  adjacency  matrix  M  in  Figure 
10,  where  the  partite  1  consists  of  nodes  {1,2,  3},  partite  2  contains  nodes  {4,  5,  6},  and  partite  3  contains 
nodes  {7,  8,  9}.  BitCLQ  is  initialized  by  setting  S  :=  0,5^  :=  {1,2,3}  and  Z\  :=  {111|111|111}. 
Since  all  the  partites  arc  of  the  same  size,  i.e.  \Zl  b\  =  3  for  all  h  e  B s,  the  one  to  branch  on  is  chosen 
arbitrarily;  assume  that  the  first  part ite  Z\  \  is  chosen  for  branching.  The  search  process  from  this  point 
restarts  3  times,  each  time  adding  one  of  the  three  nodes  in  Z\  \.  The  first  node  to  add  to  S  is  node  1, 
Z  14,1  is  then  set  to  0,  and  Z2  is  subsequently  created  by  performing  logical  AND  operation  with  Z\ 
and  the  corresponding  row  of  the  adjacency  matrix  M  as  operands: 
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Algorithm  2  BitCLQ(f) 


1 

b,  €  argminfc  {\Zuh\  \  b  e  B s} 

2 

i  :=  the  first  significant  bit  in  Zt  j,t 

3 

repeat 

4 

n  t  :=  the  i -th  node  in  bt 

5 

Zt,b,i  '■=  0 

6 

S  :=  S  U  {«,} 

7 

if  S  =  k  then 

8 

Q:=QUS 

9 

S:=S\{nt} 

10 

else 

11 

Zf  +  l,b  -=  Z /_ /;  A  M(l1t 

)  for  all  b  e  B  s 

12 

Bs  '■=  Bs  U  {bt}',  Bs 

=  Bs\{b,} 

13 

P  :=  number  of  partitions  Zt  b  with  Zt  j, 

14 

if  1 5  +  P  =  k  then 

15 

BitCLQfi  +  1) 

16 

S  :=S\{n,} 

17 

Bs  '■=  Bs\{bt}', 

Bs  '■=  Bs  U  {bt} 

18 

else 

19 

S  :=S\{nt} 

20 

CO 

cq 

II 

Co 

cq 

Bs  '■=  Bs  U  {b,} 

21 

end  if 

22 

end  if 

23 

i  :=  the  first  significant  bit  in 

Zt,bt 

24 

until  i  <  \Vht\ 

t 

=  1, 

s 

:=  {1}, 

Z2  :=  Zj  A  M{\)  =  {011|111|111}  A  {000|  1 1 1 1 0 1 1 }  =  {000|  1 1 1 101 1}, 

BS  :=  {2,3}. 

As  a  result,  the  set  Z?  of  nodes  compatible  with  the  partial  solution  S  —  {1}  contains  nodes  {4,  5,  6,  8,  9}. 
Since  none  of  the  parti tes  in  B $  is  empty,  the  partial  solution  S  is  promising  and  a  new  subproblem  is 
created.  The  objective  in  the  new  subproblem  is  to  find  a  |5s|-clique  in  Z2.  A  node  from  Z2; 3  will  be 
added  to  S  (since  | Z2>3 1  <  |Z2>2|).  The  first  node  in  Z2; 3  to  add  to  the  partial  solution  is  node  8.  The 
bit  corresponding  to  node  8  is  set  Z2p)2  :=  0,  and  we  have 

t  :=  2, 

5:={1,8}, 

Z3  :=  Z2aM(8)  =  {000|111|001}  A  { 1 1 1 1001 1000}  =  {000|001|000}, 

BS  :=  {2}. 

Again,  the  partites  in  B s  contain  at  least  1  node  (node  6)  in  Z3.  So  the  partial  solution  is  promising,  and 
a  new  subproblem  is  created.  In  the  next  step,  node  5  is  added  to  S : 

t  :=  3, 

5  :=  {1,8,6}. 

At  this  point,  since  |Sj  =  k  —  3,  i.e.,  a  /:-clique  is  found.  To  continue  the  search  for  other  A: -cliques,  the 
last  node  in  S  is  removed.  BitCLQ  searches  Z3-2  for  another  node  that  can  be  added  to  S.  Since  such  a 
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Figure  10:  A  3-partite  graph  and  its 


node  does  not  exist,  the  algorithm  backtracks:  t  :=  2,  node 
with  S  —  { 1 ,  9}  as  the  partial  solution. 

6.2  Numerical  Results 

In  order  to  illustrate  the  performance  of  the  proposed  method,  the  k -clique  enumeration  problem  for 
A' -part ite  graphs  has  been  solved  by  BitCLQ  and  FINDCLIQUE  for  randomly  generated  graph  instances 
of  several  types.  Both  algorithms  were  implemented  in  C++  and  ran  on  a  64-bit  Windows  machine 
with  3GHz  dual-core  processor  and  4GB  of  RAM.  It  is  worth  noting  that  the  original  implementation 
of  FINDCLIQUE  algorithm  by  Grunert  et  al  Grunert  et  al.  (2002)  relies  on  the  use  of  vectors  and 
links  data  types  from  the  C++  standard  template  library  (STL).  In  our  experiments,  we  observed 
that  by  replacing  the  original  data  structure  of  vectors  of  lists  with  arrays,  up  to  300%  improvement  in 
FINDCLIQUE  running  time  is  achieved  on  the  data  sets  used  in  our  case  study.  The  numerical  results 
reported  for  the  FINDCLIQUE  algorithm  are  obtained  using  this  “improved”  implementation. 

Our  numerical  experiments  involve  randomly  generated  instances  of  k -partite  graphs  of  two  types.  The 
first  set  of  instances  consists  of  two  groups:  small-size  instances  and  large-size  instances.  In  the  small- 
size  instances,  A -part ite  graphs  are  randomly  generated  with  the  number  of  partites  in  the  range  k  e 
[3,  10].  For  each  value  of  k,  the  reported  running  times  and  the  number  of  A -cliques  in  the  graph  are 
averaged  over  10  instances.  Table  9  shows  the  summary  of  the  experimental  results  for  this  first  group. 
The  columns  of  the  table  show  the  number  k  of  partites  in  the  k  -partite  graph,  the  number  in  of  nodes 
in  each  partite  of  the  graph,  the  total  number  \V\  of  nodes  in  the  graph,  the  graph’s  density  p,  and 
the  total  number  of  A -cliques  in  the  graph  (#CLQ).  The  density  parameter  p  is  used  for  generation 
of  the  graphs,  and  is  equal  to  the  probability  of  an  edge  connecting  two  nodes  from  different  partites: 
Pr  {(vi,  Vj)  e  E)  =  p. 

The  second  group  include  instances  of  larger  size  with  the  values  of  k  e  {25,  50,  75,  100}.  For  each 
value  of  k  in  this  group,  10  random  instances  of  the  A: -partite  graph  have  been  generated  and  solved  by 
FINDCLIQUE  and  BitCLQ.  Table  10  summarizes  the  results  of  the  experiments  for  this  group.  Since 
the  graphs  used  in  this  set  of  experiments  are  rather  large  and  the  list  of  all  A: -cliques  contained  in  them 
may  not  be  found  in  a  reasonable  time,  the  solution  process  has  been  terminated  after  200  seconds  and 
the  number  of  A -cliques  found  by  each  method  was  recorded.  BitCLQ  outperformed  FINDCLIQUE  in 
all  cases. 

The  second  set  of  experiments  was  conducted  to  compare  the  performance  of  BitCLQ  with  FIND¬ 
CLIQUE  on  randomly  generated  instances  of  Multidimensional  Assignment  Problem  (MAP).  As  shown 
in  Krokhmal  et  al.  (2007);  Krokhmal  and  Pardalos  (2011);  Mirghorbani  et  al.,  high-quality  solutions  for 


3  4  5  6  7  8 

0  1  1  1  0  1  1\ 

0  0  0  1  0  1  1 

0  0  0  0  0  1  1 

0  0  0  0  1  0  0 

0  0  0  0  1  0  1 

0  0  0  0  1  1  0 

0  1110  0  0 

1  0  0  1  0  0  0 

1  0  1  0  0  0  0/ 

adjacency  matrix. 


8  is  removed  from  S,  and  BitCLQ  restarts 
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Table  9:  Average  computational  time  (in  seconds)  to  find  all  the  ^-cliques  (#CLQ)  contained  in  randomly 
generated  k- partite  graphs. 


k 

m 

m 

P 

#CLQ 

FINDCLIQUE 

BitCLQ 

3 

100 

300 

0.1 

1004 

0.005 

0.002 

4 

100 

400 

0.15 

1124 

0.008 

0.002 

5 

100 

500 

0.2 

1047 

0.015 

0.003 

6 

100 

600 

0.25 

939 

0.031 

0.006 

7 

50 

350 

0.35 

192 

0.009 

0.004 

8 

50 

400 

0.4 

299 

0.021 

0.007 

9 

50 

450 

0.45 

683 

0.055 

0.021 

10 

50 

500 

0.5 

2672 

0.176 

0.071 

Table  10:  Average  number  of  /r -cliques  found  in  randomly  generated  instances  of  k -partite 
200  seconds. 

k 

m 

m 

P 

time 

FINDCLIQUE 

BitCLQ 

25 

40 

1000 

0.8 

200 

13,556,733 

23,516,581 

50 

30 

1500 

0.9 

200 

800,369 

1,032,111 

75 

30 

2250 

0.95 

200 

557,042,389 

735,722,241 

100 

30 

3000 

0.95 

200 

348,416 

365,799 

randomized  MAPs  can  be  obtained  as  n -cliques  in  n  -partite  graphs  that  are  constructed  in  a  special  way 
from  the  problem’s  data  (in  this  case,  n  denotes  the  number  of  elements  per  dimension  in  a  d  -dimensional 
MAP).  For  MAPs  with  random  iid  costs,  the  resulting  n -partite  graph  can  be  viewed  as  randomly  gen¬ 
erated  with  a  certain  density.  The  corresponding  results  arc  reported  in  Table  1 1 ,  where  n  denotes  the 
number  of  partitions  in  the  graphs,  and  cl  is  the  number  of  dimensions  d  in  the  MAP. 
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